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Hybrid  Finite  Element/Waveguide 
Mode  Analysis  of  Passive  RF  Devices 

1.  INTRODUCTION  AND  BACKGROUND 

A  numerical  solution  for  the  time-harmonic  electric  fields  inside  a  cavity  fed  by  two  waveguides 
at  opposite  ends  has  been  implemented  and  validated.  The  generic  problem  is  that  of  a  two-port  RF 
(radio  frequency)  circuit.  Specific  test  results  were  obtained  for  a  variety  of  problems,  including  wave¬ 
guide  discontinuities  and  microwave  circuits.  The  emphasis  of  this  report  is  on  the  analysis  technique, 
which  uses  the  finite  element  method  in  conjunction  with  waveguide  modes  in  a  hybrid  scheme.  The 
results  of  this  work  are  both  a  demonstration  of  the  capabilities  that  may  be  obtained  using  the  finite 
element  method  and  a  computer  code  that  is  accurate  and  versatile. 

This  is  an  interim  report  in  a  project  dealing  with  the  analysis  of  antennas  using  the  finite  element 
method  (FEM).  The  results  described  here  demonstrate  some  of  the  key  building  blocks,  especially  the 
coupling  of  the  FEM  solution  to  higher-order  waveguide  modes.  This  will  lead  to  an  accurate  feed 
description  for  antennas  that  are  driven  by  waveguides.  There  are  other  methods  available  for  waveguide 
discontinuity  and  RF  circuit  analysis,  but  the  present  approach  is  partly  justified  by  its  intended  applica- 
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tion.  Nonetheless,  the  novel  approach  of  using  higher-order  waveguide  modes  in  conjunction  with  the 
finite  element  method  has,  for  many  problems,  a  decided  advantage  over  previous  solutions  that  used  only 
the  lowest  order  waveguide  mixie  and  extended  the  finite  element  region  into  the  waveguide  far  enough 
to  make  evanescent  inodes  negligible. 

This  report  will  discuss  the  general  cavity/waveguide  problem  and  its  formulation  as  a  variational 
problem  using  weighted  residuals  and  Galerkin's  method  (Sections  2  and  3).  These  two  chapters  are 
important  to  understanding  how  the  physical  problem  is  cast  as  a  mathematical  boundary  value  problem. 
Section  4  will  show  the  details  of  the  derivation  of  a  matrix  equation  using  vector  edge-based  finite 
elements,  as  well  as  the  methods  used  to  solve  the  equation  and  obtain  S  parameter  information  from  the 
solution.  That  chapter  is  important  to  the  reader  who  wishes  to  understand  the  details  of  the  finite 
element  and  waveguide  mode  solutions.  Section  5  will  discuss  the  structure  of  the  analysis  code,  includ¬ 
ing  its  interface  to  commercial  software  used  for  geometry  definition  and  mesh  generation.  That  informa¬ 
tion  will  be  useful  mainly  to  the  reader  interested  in  using  or  modifying  the  code,  or  in  developing  a 
similar  code.  Last,  Section  6  will  show  the  results  of  test  cases  that  validate  the  code’s  accuracy  and 
demonstrate  its  usefulness. 

2.  PROBLEM  DESCRIPTION  AND  SOLUTION  APPROACH 

The  general  problem  to  be  solved  is  illustrated  in  Figure  1 .  The  cavity  region,  denoted  0,  has 
perfectly  conducting  walls  except  for  apertures  at  z=0  and  z=d.  The  apertures  need  not  extend  over  the 
entire  waveguide  cross  sections.  The  boundary  of  0  is  denoted  P,  and  the  two  apertures’  open  areas  are 
and  Pg.  The  two  waveguides  need  not  be  the  same  size  or  even  the  same  type.  Derivations  and 
results  will  be  presented  for  coaxial,  circular  and  rectangular  waveguides.  The  apertures  are  in  the 
waveguides’  end  walls. 

The  cavity  interior  may  be  regarded  as  an  arbitrary  two-port,  passive,  linear  device.  It  may 
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Waveguide  A 


Cavity 


Waveguide  B 


2  =  0  z  =  d 

Figure  1.  Generic  Waveguide/Cavity  Problem 

contain  any  mixture  of  isotropic  dielectric  and  ferrite  materials,  as  well  as  perfectly  conducting  obstacles. 
Those  features  are  adequate  to  describe  a  broad  range  of  practical  RF  devices.  Some  examples  are  shown 
in  Figure  2.  In  the  first  (Figure  2a)  the  inlet  waveguide  is  a  coaxial  line  whose  center  conductor  extends 
into  the  cavity  region,  a  short  section  of  rectangular  waveguide.  The  outlet  waveguide  is  the  same  size 
rectangular  waveguide.  In  the  second  example  (Figure  2b)  the  cavity  is  a  rectangular  conducting  box 
containing  a  dielectric  slab  that  supports  a  microstrip  conductor  (a  meander  line).  The  inlet  and  outlet 
ports  are  both  coaxial  waveguides. 

An  overview  of  the  solution  procedure  is  as  follows:  (» /  Represent  the  waveguide  fields  in  terms 
of  orthonormal  mode  functions;  (2)  Represent  the  field  inside  0  by  a  variational  functional;  (3)  Match 
the  transverse  magnetic  field  across  the  apertures  by  incorporating  the  waveguide  mode  fields  in  the 
functional’s  boundary  term;  (4)  Generate  a  linear  system  by  expanding  the  field  in  term;  of  known 
functions;  and  (5)  compute  the  devices’  S  parameters  from  the  solution  of  the  linear  system  (matrix 
equation).  The  following  two  sections  provide  the  details  of  steps  1  -  3  and  4  &  5,  respectively. 


3.  VARIATIONAL  STATEMENT  AND  BOUNDARY  CONDITIONS 

3.1.  Variational  Statement 


The  finite  element  method  is  usually  applied  to  a  functional,  an  integral  that  contains  the  unknov  n 
field  as  part  of  the  integrand  (Ij.  (A  functional  is  a  mapping  of  a  function  into  a  number.  In  other 
words,  the  result  of  the  integration  is  a  single  number.  This  is  in  contrast  to  an  integral  equation,  which 
al.so  contains  an  unknown  function  in  the  integrand,  but  which  results  not  in  a  number,  but  in  another 
function.)  The  functional  is  also  a  variational  statement,  meaning  that  the  correct  solution  for  the 
unknown  field  is  the  function  that  minimizes  the  integral.  In  this  case,  the  variational  statement  is  a 
form  of  the  time-harmonic  wave  equation  for  the  electric  field  inside  fi.  The  source-free  form  is  appro¬ 
priate  to  this  problem  because  all  sources  are  assumed  to  be  located  at  z<0  in  waveguide  A.  The  weak 
form  is  the  functional 


F(E)  =  I 


()  L 


Mr 


.  V  X  W *  •  V  xE  -k-Qt^W*  •  E 


dv 


-  J^O^O  I  =  0  (1 


where  E  is  the  unknown  electric  field,  W  is  a  trial  function,  H  is  the  magnetic  field,  is  the  free  space 
wavenumber,  t/q  is  the  free  space  impedance  and  fi  is  the  outward  normal  to  F.  This  is  applicable  to 
linear,  isotropic,  inhtimogeneous  materials.  Appendix  A  discu.sses  how  this  is  derived  and  why  this  form 
is  more  appropriate  than  others  that  appear  in  the  electromagnetics  literature.  This  functional  is  a 
statement  of  conservation  of  energy,  with  the  volume  integral  representing  stored  energy  and  dissipation 
in  materials,  and  the  surface  integral  representing  power  flow  into  and  out  of  the  volume.  For  that 
reason,  the  perfectly  conducting  portions  of  F  do  not  contribute  to  the  surface  integral,  and  it  reduces  to 
an  integration  over  F^  and  F^,  where  it  will  be  used  to  enforce  continuity  of  tangential  H. 
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3.2.  Waveguide  Field  Representation 

The  following  derivation  is  appropriate  for  arbitrary  waveguide  modes.  It  is  valid  for  any  type 
of  waveguide  for  which  orthonormal  mode  functions  are  known:  rectangular,  circular,  coaxial  and 
elliptical.  Appendix  B  gives  the  specific  forms  for  coaxial  and  circular  waveguides.  Appendix  C 
discusses  rectangular  waveguides. 

In  waveguide  A,  the  electric  field  transverse  to  z  may  be  written  in  terms  of  orthonormal  mode 
functions  [2), [3]; 

1=0 

The  index  i  includes  TE  and  TM  modes  (and  TEM  if  applicable)  with  /=0  corresponding  to  the  dominant 
mode.  The  first  term  represents  the  incident  field,  assumed  to  be  of  unit  magnitude  and  in  the  dominant 
mode  only.  The  unknown  coefficients  C,  are  the  excitation  coefficients  for  the  modes  reflected  by  the 
device.  Cg  is  particularly  significant:  Cg-l  is  S||,  the  reflection  coefficient  seen  by  an  observer  in  wave¬ 
guide  A.  A  similar  expression  holds  for  waveguide  B,  except  that  there  is  no  incident  field  term,  and 
the  coefficient  for  the  dominant  mode  will  be  proportional  to  S21,  the  transmission  coefficient.  The  mode 
functions  are  entirely  transverse  to  the  waveguide  axis  (the  z  axis).  The  propagation  constants  are  related 
to  the  mode  cutoff  wavenumbers  by 


7,  = 


(3) 


which  is  positive  imaginary  for  propagating  modes  and  positive  real  for  evanescent  modes. 
The  transverse  magnetic  field  in  waveguide  A  is 
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(4) 


-  £  C,Y,e^‘Ut^g,) 
1=0 


where  K,  is  the  modal  admittance; 


-  ‘  (TE) 
io’tOt'r 

Y-  =  ^  (5) 

yino 

i;*'  (TEM) 


The  fields  in  waveg'oide  B  are  essentially  identical,  except  that:  (a)  there  is  no  incident  field  tern;  (b)  the 


fields  are  propagating  in  the  +z  direction;  and  (c)  its  mode  functions,  modal  admittances,  and  propagation 


constants  need  not  be  the  same  as  waveguide  A. 


3.3.  Transverse  Field  Continuity  at  the  Apertures 


At  2=0  the  outward  normal  to  fl  is  -z,  so  from  Eq.(4); 


=  roSo  (6) 

z=0  ,,0 


The  Cj’s  are  found  by  taking  the  inner  product  of  Eq.(2)  with  every  over  the  waveguide  cross  section 
and  using  the  orthogonality  property,  with  the  result: 


The  limits  of  integration  above  are  over  the  waveguide  end  wall  aperture  instead  of  its  cross  section 
because  the  transverse  electric  field  is  zero  on  the  conducting  surfaces  at  z  =  0. 

Now  let  (F^-Finc)  represent  that  part  of  the  boundary  integral  in  Eq.(l)  over  and  substitute 
Eqs.(6)  and  (7): 
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-jkovol  -2  Ko  J  VV  *  •  gods  -  E  I  ^  •  k,ds  I 


(9) 


The  first  term  inside  the  braces  in  Eq.{9)  is  the  negative  of  the  incident  field  term, 

In  the  same  manner,  the  boundary  functional  for  waveguide  B  can  be  written  in  terms  of  its  mode 
functions  /i,  and  modal  admittances  f/: 


00 

=  yl^oVoT,  K  f  •  hids  I  if-  lids 


Let  F/  denote  the  volume  integral  in  Eq.(l),  which  is  the  part  of  the  functional  dealing  with  the 
interior  fields.  Then  the  functional  equation  may  be  written  symbolically  as 

(ID 

When  the  unknown  and  trial  fields  are  expanded  in  series  of  known  functions  with  unknown  coefficients, 
Eq.(l  1)  will  become  a  matrix  equation; 

=  [5^  +  5'^ +5®]£ 

where  £*"‘^  is  a  column  vector  representing  the  incident  field,  £  is  a  column  vector  representing  the 
unknown  field  at  locations  inside  the  volume,  and  Sf,  and  5®  are  matrices.  The  following  section 
gives  the  details  of  the  discretization  process  by  which  the  continuous  functional  is  mapped  into  the  linear 
system  to  give  expressions  for  the  entries  in  £‘"‘',  5^  5^  and  5®  in  terms  of  known  quantities. 
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4.  DISCRETIZATION 


The  finite  element  method  generates  a  linear  system  of  equations  from  a  linear  functional.  The 
alternative  used  here  is  Galerkin’s  method,  for  reasons  explained  in  Appendix  A.  It  is  quite  similar  to 
the  method  of  moments,  the  distinction  being  that  MOM  applies  weighted  residuals  to  an  integral  equation 
rather  than  to  a  functional  (Reference  4,  p.  18). 

The  first  three  parts  of  this  section  discuss  the  expansion  functions  that  were  chosen  to  represent 
the  field  and  trial  functions:  the  relative  merits  of  scalar  and  vector  expansion  functions;  the  structure 
of  the  scalar  functions  and  their  representation  in  homogeneous  coordinates;  and  the  construction  of  the 
vector  expansion  functions  from  the  scalar  functions.  The  fourth  part  discusses  the  weighted  residuals 
procedure  leading  to  the  matrix  equation.  The  fifth  and  sixth  parts  derive  explicit  expressions  for  the 
matrix  elements  and  the  incident  field  terms.  The  last  part  discusses  methods  for  solving  the  matrix. 

4.1.  Scalar  vs  Vector  Expansion  Functions 

Figure  3  illustrates  a  gridding  of  a  volume  region  into  small  elements.  These  elements  or  cells 
are  tetrahedra,  each  of  which  has  four  faces  and  four  vertices,  or  nodes.  This  figure  illustrates  an 
important  feature  of  the  finite  element  method:  the  grid  density  may  be  varied  within  the  problem  space 
so  that  it  is  coarse  in  regions  where  the  fields  are  expected  to  be  well-behaved;  and  fine  in  regions  where 
the  geometry  is  more  complex  or  where  fields  may  have  singularities. 

The  electric  field  inside  the  region  shown  is  to  be  expanded  in  terms  of  some  form  of  known 
functions  defined  over  individual  cells  with  unknown  coefficients  representing  the  field  strength  and 
direction.  The  decision  of  whether  to  use  scalar  functions  with  vector  coefficients  or  vector  functions 
with  scalar  coefficients,  as  contrasted  in  the  following  two  equations,  is  much  more  than  a  matter  of 
preference. 
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Figure  3.  Subdivision  of  a  Volume  Region  (one  quadrant  of  a  disk)  Into  Tetrahedra; 
(a)  interior  edges  visible;  (b)  interior  edges  hidden 


M 


f=i 


(13) 


N 

£■  =  52 

5  =  1 

Traditional  finite  element  calculations  involve  solving  for  the  unknown  function  at  the  grid  nodes. 
Hence,  Eq.(13)  seems  the  most  natural  expansion  for  the  electric  field.  M  is  the  total  number  of  nodes 
in  the  volume.  The  function  is  a  subdomain  function,  meaning  that  it  is  only  defined  over  a  small 
portion  of  the  volume,  specifically,  those  cells  adjacent  to  node  s.  In  contrast,  the  index  s  in  Eq.(14)  is 
over  the  grid  edges  (the  lines  joining  nodes). 

The  node-based  expansion  in  Eq.(13)  has  three  important  disadvantages: 

(a)  It  is  difficult  to  enforce  the  boundary  condition  flxE=0  at  nodes  on  perfect  conductors 
(mainly  in  an  algorithmic  sense).  Worse  yet,  if  a  node  is  on  a  conducting  edge  or  tip,  it  is  not 
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possible  to  enforce  the  condition  because  A  is  not  defined.  The  only  remedy  is  to  refine  the  mesh 
near  edges  so  that  their  erroneous  contributions  are  less  significant. 

(b)  The  expansion  of  Eq.(13)  does  not  meet  the  divergence  condition  V-(e£)  =  0.  It  has  been 
presumed  that  this  could  be  circumvented  by  using  a  penalty  function,  but  Boyse  et  al.  point  out 
that  penalty  methods  are  only  justified  for  positive  definite  functionals  [5],  whereas  Eq.(l)  is 
indefinite. 

(c)  Suppose  that  the  node  s  is  on  an  interface  between  dielectrics.  Then  the  component  of  £ 
normal  to  the  interface  must  be  discontinuous.  Yet  Eq.(13)  implies  that  all  components  of  E  are 
continuous  at  s.  Hence  the  scalar  functions  over-prescribe  the  continuity  at  dielectric  interfaces 
(6).  The  remedy  is  to  decrease  the  mesh  cell  size  around  dielectric  boundaries  so  that  the  error 
contributed  by  the  boundary  nodes  is  less  significant. 

All  three  of  these  disadvantages  are  sources  of  inaccuracy  and  all  three  incur  a  computational  cost:  (a) 
and  (c)  in  terms  of  matrix  size:  and  (b)  in  terms  of  matrix  element  calculation  time. 

Several  authors  have  reported  success  using  vector  expansion  functions  as  in  Eq.(14)  in  cases 
where  the  scalar  fonctions  failed  |6],(71.  The  specific  form  of  the  function  will  be  described  in  detail 
in  Section  4.3.  Its  advantages  are: 

(a)  It  meets  the  divergence  condition. 

(b)  It  enforces  only  tangential  field  continuity  across  boundaries;  and 

(c)  It  allows  a  simple  yet  effective  method  for  implementing  conductor  boundary  conditions, 
even  at  edges. 

The  first  two  are  vital  in  terms  of  using  the  finite  element  method.  These  two  qualities  assure  that  the 
expansion  and  trial  functions  are  admissible.  The  vector  element’s  disadvantages  are: 

(a)  It  is  not  compatible  with  most  existing  preprocessing  and  postprocessing  software  (usually 
CAD  packages  oriented  to  finite  element  computations  for  mechanical  engineering),  so  additional 
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algorithms  and  computations  are  needed  to  transform  the  input  geometry  and  output  results;  and 
(b)  There  are  generally  4-5  times  as  many  edges  as  nodes  in  a  tetrahedron  mesh,  and  so  the 
edge-based  matrix  will  be  4/3  -  5/3  larger.  This  is  partly  ameliorated  by  the  fact  that  there  is  less 
need  for  fine  sampling  at  material  interfaces  and  conductor  edges.  Furthermore,  the  connectivity 
of  edges  is  less  than  for  nodes,  so  the  edge-based  matrix  is  considerably  more  sparse,  and  thus 
the  total  number  of  matrix  entries  may  actually  be  less.  The  overall  matrix  size  is  a  greater 
concern  when  using  LU  decomposition,  but  the  total  number  of  nonzero  entries  is  the  most 
important  for  an  iterative  solver  such  as  the  conjugate  gradient  method. 

In  summary,  the  edge-based  vector  functions  were  adopted  because  their  disadvantages  affect  algorithm 
complexity,  while  their  advantages  in  theoretical  and  mathematical  correctness  provided  a  path  with 
considerably  less  risk. 

4.2.  Linear  Expansion  Functions  on  Triangles  and  Tetrahedra 

Since  the  vector  expansion  functions  ^  are  defined  in  terms  of  the  more  conventional  scalar 
functions  <^,  the  latter  are  described  here  in  some  detail.  Each  is  defined  over  the  cells  (tetrahedra) 
adjacent  to  node  s.  Its  value  is  1  at  node  s,  and  zero  at  the  surrounding  nodes  and  opposing  faces.  It 
may  in  general  have  a  polynomial  variation  over  the  cells  in  which  it  is  defined,  but  the  present  work  uses 
linear  functions  only. 

A  two-dimensional  linear  function  is  shown  in  Figure  4.  It  is  defined  over  the  triangles  adjacent 
to  a  single  node,  where  its  value  is  1.  This  function  may  be  regarded  as  being  assembled  from  "linear 
finite  elements"  that  are  defined  on  individual  triangles.  Each  triangle  supports  three  finite  elements,  each 
of  which  is  1  at  a  single  node  and  zero  along  the  opposite  edge.  Figure  5  illustrates  the  three-dimensional 
equivalent:  The  node  is  at  the  center  of  a  group  of  8  tetrahedra  that  form  a  diamond  shape.  One  tetrahe¬ 
dron  is  shown  expanded  with  dashed  lines  indicating  surfaces  where  the  function  has  constant  value.  The 
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Figure  4.  A  Two  Dimensional  Linear  Expansion  Function  (defined  over  triangular 
cells  adjacent  to  a  node  and  zero  on  all  other  triangles) 


Figure  5.  A  Three  Dimensional  Linear  Expansion  Function:  (left)  eight  cells 
surrounding  a  node;  (right)  surfaces  of  constant  function  value  in  one  cell 
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three-dimensional  is  constructed  from  functions  that  are  defined  over  individual  tetrahedra.  For 
example,  within  the  tetrahedron  with  index  k,  there  are  four  linear  finite  elements  one  for  each  node 
/■=  1, 2, 3,4.  The  index  i  is  a  local  node  index  while  s  is  a  global  node  index.  An  expression  for  <\>^  is 

EE y;'*’ «(*■;'•*)  <'« 

i  =  l  1  =  1 

The  notation  5(s:i.k)  is  a  Kronecker  delta  that  equals  1  if  the  global  node  s  is  the  same  as  the  local  node 
i  within  cell  k,  and  equals  zero  otherwise.  K  is  the  total  number  of  cells  in  the  mesh.  The  notation  in 
Eq.(15)  is  merely  a  formal  representation  for  the  assembly  procedure  that  allows  matrix  element  calcula¬ 
tions  to  be  performed  one  cell  at  a  time,  which  is  considerably  more  efficient  than  performing  them  one 
node  at  a  time. 

The  integrations  leading  to  matrix  elements  are  products  of  these  linear  functions  and  their  partial 
derivatives.  Those  integrations  may  be  performed  easily  in  closed  form  using  homogeneous  coordinates 
(also  known  as  simplex  coordinates)  [Reference  8,  pp.266-274J.  The  coordinate  transformation  is 
defined  locally  within  each  cell.  In  the  tetrahedron  there  are  four  coordinates  r,,  h  ^4- 
coordinate  /,  of  a  point  anywhere  within  the  cell  is  the  normalized  distance  to  node  /  from  the  opposing 
face  (normalized  to  the  height  of  the  cell  in  that  direction).  The  coordinate  transform  is  given  in  terms 
of  a  matrix  T. 


where  V  is  the  cell  volume.  The  entries  in  T  are  the  16  cofactors  of  the  following  matrix  made  up  of  the 
cell  vertex  coordinates: 
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1  Jr,  y,  z, 

I  Xy  y-,  Z-, 

U  =  -  -  - 

J  -*3  >'3  23 

1  J:4  ^4  24 

The  two  reasons  these  coordinates  are  so  convenient  are:  (a)  the  linear  expansion  functions  are 

simply 

and  (b),  the  limits  of  integration  are  simplified,  leading  to  the  following  result  for  the  inner  product  of 
any  two  linear  expansion  functions  over  a  single  cell: 


I  I  \fifjdxdydz 

cell  k 


6V 


I  r//,  J  dh  J  t,tjdt^ 
0  0  0 


(19) 


This  simplicity,  one  of  the  attractive  features  of  FEM,  applies  equally  to  the  vector  expansion  functions. 
4.3.  Vector  Edge-Based  Expansion  Functions 

The  vector  expansion  function  is  defined  relative  to  the  edge  joining  nodes  i  and  j  [7]: 

fij- ™ 
where  and  fj  are  the  linear  node-based  functions  and  L^j  is  the  edge  length.  The  scaling  by  L^j  ensures 
that  the  component  of  tangential  to  the  edge  is  a  unit  vector.  Figure  6  illustrates  the  two  dimensional 
version  of  this  function.  Observe  that  Vf^  is  a  vector  in  the  direction  of  steepest  ascent  of  so  it  is 
directed  inward  from  the  face  opposite  node  /.  This  choice  of  function  has  the  important  properties  that: 
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^  2 

/> 

'■  ■/ 

2  FUNCTION  DEFINED 


Figure  6.  A  Two-Dimensional  Vector  Expansion  Function:  (a)  direction 
(arrow  length  indicates  function  magnitude);  (b)  constituent  linear 

finite  elements 


(a)  V  •  so  the  divergence  condition  is  automatically  satisfied;  and  (b)  the  component  along  the  edge 
i-j  has  constant  unit  magnitude,  ensuring  tangential  field  continuity  across  element  boundaries.  One  can 
represent  the  global  expansion  functions  as  being  assembled  from  functions  defined  over  individual  cells; 

"  E  E  E  ^fj  v,(^))  (21) 

*=1  i=I  y=l 

where  8{s;iJ,{k))  is  a  Kronecker  delta  equal  to  1  when  global  edge  s  is  the  same  as  the  edge  connecting 
local  nodes  i  and  j  within  cell  k.  This  completes  the  definition  of  the  vector  edge-based  functions.  Their 
use  in  the  discretization  of  the  functional  is  discussed  next. 


4.4.  Discretization  Using  Galerkin’s  Method 


Galerkin’s  method  is  the  specialization  of  weighted  residuals  in  which  the  trial  functions  are  the 
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same  as  the  expansion  functions.  Let  be  the  weighting  function  for  edge  r.  Then  substituting  tV»  = 
for  each  r=  into  (1)  will  generate  N  simultaneous  equations 


Fr 


Mr 


dv 


-  jkQi]o^^r  =  Q,  r=l,2 . N 


(22) 


(The  subscript  r  on  /x  and  e  still  denotes  relative  permeability  and  permittivity.)  Substituting  the  expan¬ 
sion  for  E  from  Eq.(l4)  into  the  two  terms  of  the  volume  integral: 

N  r  1 


/.r 


p/;  t./2 


(23) 


5=1 


f  _L  V  x^  •  V  dv 
1 
0 


(24) 


(25) 

Q 

To  simplify  the  surface  integral  terms,  let  denote  the  inner  product  of  the  r’th  expansion  function  and 
the  i’th  waveguide  mode  function: 


(26) 


<,  =  j 

then  from  the  First  term  of  Eq.(9) 

4"  =  VWo^rO 


(27) 


(28) 
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and  from  the  second  term, 


•yio'IoE  y.*r,*,. 

1=0 


(29) 


and  similarly,  from  Eq.(lO) 

00 

s’  -J^oloE 

i=0 


(30) 


Collecting  the  terms  from  Eqs.(24)-(30)  gives  the  matrix  equation  (12).  The  next  two  sections  give  the 
details  for  performing  the  integrations  in  Eqs.(24)-(27). 


4.5.  Volume  Integral  Computations 


The  volume  integral  computations  are  carried  out  by  visiting  each  cell  once  and  adding  its 
contribution  to  S„  for  each  combination  of  edges  shared  by  the  cell,  excluding  those  edges  that  are  on 
perfect  conductors.  Let  ij  and  m,n  be  the  local  node  indices  defining  edges  r  and  s,  respectively, 
1  <ij,m,n<4,  i^j,  m^n.  After  carrying  out  the  coordinate  transformation  local  to  cell  k,  giving  the 
matrix  T: 


(31) 


(32) 


Using  the  identity  Vx(aVb)  =  aVxVb  +  VaxVb  and  noting  that  the  second  derivatives  of  the  linear 
functions  /are  all  zero. 


=  2L,Vf,xVfj 
Vx(f,  =  2L,V/^xV/„ 


(33) 
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The  terms  in  Eq.(33)  are  all  constants  and  may  be  taken  outside  the  integral  in  Eq.(24) 


^  [  (7;37-4  - 

*  ^^i2^3"^3^2)(^m2^n3~^m3^n2)] 


(34) 


The  second  volume  integral’s  contribution  from  cell  k  is 


Sll  =  -klt^L^L^  J  •  V/„ .  V/„ 

‘ -yi/„v/,  •  v/„  .  v/„]rfv 


(35) 


After  breaking  this  integral  into  four  parts,  the  gradient  terms  are  constants  that  may  be  brought  outside 
the  integrals,  which  may  then  be  evaluated  using  Eq.(19).  The  result  is 


‘"’''‘'‘'Elo  7;, 


720  V 


/=2 


(36) 


In  Eqs.(34)  and  (36),  the  volume  integral  has  been  evaluated  in  closed  form  and  written  as  an  algebraic 
expression  in  terms  of  the  geometry  of  a  cell  and  the  constitutive  parameters  contained  within  it.  The 
following  section  derives  similar  expressions  for  the  matrices  representing  the  waveguide  apertures. 


4.6.  Boundary  Integral  Computations 


The  only  difficult  computation  involved  in  the  boundary  integrals  Eqs.(26)-(28)  are  the  inner 
products  When  the  functions  are  circular  or  coaxial  mode  functions  the  integrals  may  not  be 
evaluated  in  closed  form,  so  they  must  be  evaluated  numerically. 

As  was  the  case  with  the  volume  integrals,  the  contributions  to  the  surface  integrals  are  carried 
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out  one  cell  at  a  time.  The  two-dimensional  simplex  coordinates  are  used  because  the  2D  analog  of  the 
vector  expansion  function  evaluated  over  a  triangular  tetrahedron  face  is  the  same  as  the  3D  function 
evaluated  in  that  face.  Hence  the  matrix  T gives  the  coordinate  transformation  local  to  a  triangle  and  the 
scalar  expansion  functions  are 


f.  - 


(37) 


where  A  is  the  triangle’s  area.  The  vector  functions  are  still  as  defined  by  Eq.(20),  except  that  the  index 
i  may  only  take  on  values  1,2  and  3. 

The  contribution  of  a  triangle  to  will  be  approximated  using  Gaussian  quadrature  integration. 
Quadrature  formulas  only  apply  strictly  to  one  dimension,  but  are  easily  extended  to  two  dimensional 
integration  over  rectangular  areas.  To  use  these,  the  triangle’s  geometry  is  transformed  to  a  unit  square, 
using  an  approach  suggested  by  Stroud  Sic  Secrest  (9|.  Note  that  the  transformation  to  simplex  coordi¬ 
nates  mapped  the  arbitrary  triangle  into  one  whose  vertex  coordinates  are  (0,0),  (0,1)  and  (1,0)  in 
coordinates.  The  second  transformation  is  given  by 


Ui  =  /. 


Un  = 


(I  -/,) 


.  1 


(38) 


1  ,  =  I 


The  Jacobian  of  this  transformation  is  (l-t|)  ',  so  that  a  typical  integral  term  transforms  as  follows: 


;  >-'i 


I  I 


h  =  I  I  (\  -Un)- g(Ui,U2)dU2 


(39) 
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Let  and  denote  the  une-dimensiunal  quadrature  sample  points  along  u^  and  U2,  respectively,  with 
Wj  and  the  corresponding  weights.  Then  the  integral  is  approximated  by  the  sum 

Q  Q 

k - 1  m=  I 

where  Q  is  the  order  of  the  quadrature  formula.  See,  for  example,  Reference  9  or  Reference  10,  p.887 
for  tables  of  weights  and  quadrature  points. 

Finally,  once  has  been  found  for  each  edge  s  in  for  mode  /,  then  that  mode’s  contribution 
to  may  be  computed  using  Eq.(27).  The  computation  of  is  identical.  For  rectangular  waveguides 
the  integrations  have  been  carried  out  analytically  and  appropriate  formulas  are  given  in  Appendix  C. 

The  incident  field  term  computation  is  also  the  same,  except  that  it  involves  only  the  dominant 
mode  term,  It  is  also  convenient  to  save  and  for  use  later  in  finding  the  reflection  and  trans¬ 
mission  coefficients.  Solving  the  matrix  gives  the  edge  field  coefficients  e^,  from  which  the  reflection 
coefficient  in  waveguide  A  and  transmission  coefficient  in  waveguide  B  are: 

-E -I 

J=l 

J=I 

or  Sjj  and  S^i,  respectively. 

This  completes  the  derivations  for  the  expressions  needed  to  perform  the  matrix  fill  computations. 
The  next  section  will  discuss  how  they  are  implemented  algorithmically  in  the  computer  code. 

4.7.  Matrix  Structure  and  Solution  Methods 

A  few  notes  about  the  structure  of  the  matrix  S  are  important:  (a)  Sj.^  and  Sj.^  are  zero  unless 
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edges  r  and  s  are  shared  by  one  or  more  cells,  so  they  are  sparse  matrices;  (b)  and  are  zero 
unless  edge  r  is  in  or  Fg  respectively,  so  that  the  entries  in  and  are  zero  except  when  edges  r 
and  s  are  both  in  the  same  aperture;  and  (c)  the  incident  field  elements  are  zero  except  for  those  edges 
r  in  r^.  Suppose  that  the  edges  in  the  mesh  are  indexed  sequentially  starting  with  those  in  aperture  A, 
continuing  by  increasing  z  coordinate,  and  ending  with  those  in  aperture  B.  Then:  the  matrix  S'  will  be 
banded;  5^  will  be  zero  except  for  a  dense  submatrix  in  the  upper  left  corner  of  size  (N^  is  the 

number  of  edges  in  aperture  A);  S*  will  be  zero  except  for  a  dense  submatrix  in  the  lower  right  corner 
of  size  Ng^Ng  (Ng  is  the  number  of  edges  in  aperture  B);  and  only  the  first  elements  of  are 
nonzero. 

Figure  7  shows  the  structure  of  a  typical  matrix,  with  black  points  indicating  nonzero  entries. 
Typical  of  most  problems,  the  matrix  is  banded,  but  its  bandwidth  is  quite  large.  Consequently,  there 
is  little  to  be  gained  by  using  a  band  storage  or  skyline  storage  format,  since  the  number  of  zeros  inside 
the  band  and/or  profile  is  still  quite  large.  With  a  “direct"  solver  such  as  LU  decomposition,  the  amount 
of  storage  needed  may  be  the  limiting  factor  rather  than  the  time  needed  to  perform  the  decomposition. 
Even  though  most  of  the  entries  in  S  are  zero,  the  same  is  not  true  of  its  L  (lower)  and  U  (upper)  factors 
since  there  is  a  considerable  amount  of  "fill  in"  during  the  factorization  (1 IJ.  The  problem  that  gave 
rise  to  Figure  7  had  964  unknowns.  The  matrix  had  40,462  entries  (4.4  percent  full)  while  the  LU 
factors  had  348,180  entries  (37.5  percent  full).  Heath  et  al.  review  the  necessary  steps  required  to 
efficiently  use  LU  decomposition  on  this  kind  of  matrix  112);  (a)  execute  a  reordering  algorithm  that 
swaps  rows  and  columns  in  an  attempt  to  minimize  the  fill  in  during  the  subsequent  factorization;  (b) 
perform  symbolic  LU  factorization  to  determine  where  the  fill  in  will  occur  so  that  memory  may  be 
allocated  for  the  new  matrix  elements;  and  (c)  perform  the  factorization.  But  the  reordering  algorithms 
only  lessen  the  amount  of  fill-in;  they  do  not  guarantee  that  it  will  be  less  than  any  percentage  of  the 
original  matrix  size.  Hence,  it  is  not  generally  possible  to  predict  how  much  storage  will  be  required  for 
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Figure  7.  Sample  Matrix  Structure  (dark  points  are  nonzero  elements) 


a  given  problem  aside  from  the  upper  bound  of  words. 

To  be  able  to  handle  larger  problems,  we  chose  to  use  the  conjugate  gradient  algorithm.  It  uses 
a  minimum  of  storage  because  it  always  works  with  the  original  matrix  and  does  not  generate  any  factors, 
and  therefore  there  is  no  fill-in.  The  essential  algorithm  is  given  by  Sarkar  and  Arvas  ( 13).  The  proce¬ 
dure  to  solve  AX  =  B  for  the  unknown  vector  X  begins  with  an  initial  guess  Xq,  giving  an  initial  residual 
error  vector  Rq  =  B  -  AXq.  The  initial  search  direction  and  gradient  vectors  Pq  and  Gq  are 

Po  =  Go  = 

(The  superscript  H  denotes  Hermitian,  that  is,  conjugate  transpose.)  The  iteration  consists  of  the  follow¬ 
ing  steps  using  ||  •  ||  to  denote  the  Lo  norm: 


(44) 
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(45) 

^♦1  = 

(46) 

Gi*\  =  Pi*\ 

(47) 

‘  \\G,P 

(48) 

’i*l  ^  *^iPi 

(49) 

A  typical  convergence  test  is  to  examine  the  ratio  e  =  ( ||  /?,  1|  “  /  ||  Rq  ||  at  each  iteration.  A  value  of 
10"^  or  10'^  is  usually  adequate  to  ensure  convergence  of  the  S  parameters.  This  is  demonstrated  in 
Figure  8,  which  shows  the  value  of  both  IS,,  |  and  e  for  a  test  problem  consisting  of  a  microstrip  circuit 
with  coaxial  inlet  and  outlet  waveguides.  The  reflection  coefficient  magnitude  has  converged  to  one  part 
in  1000  when  the  residual  norm  is  0.001  of  its  inital  value  (using  Xo=0  as  the  initial  guess).  It  is  also 
important  to  note  that  the  residual  norm  decreases  monotonically,  which  is  another  important  attribute 
for  a  good  convergence  criterion.  The  same  is  not  true  of  the  S  parameters  themselves  which  sometimes 
oscillate  about  their  final  values.  For  this  particular  problem  the  convergence  was  relatively  slow,  taking 
about  6N  iterations,  where  N  was  the  number  of  unknowns.  N  iterations  is  more  typical,  and  this  could 
probably  be  improved  substantially  using  preconditioning. 

5.  COMPUTER  CODE  DESIGN 

The  computer  program  ’.mplementing  the  solution  procedure  is  intended  not  only  to  demonstrate 
the  validity  and  potential  usefulness  of  the  method,  but  also  to  serve  as  a  design  tool  for  a  certain  class 
of  RF  devices.  As  such,  it  must  be  effective,  reliable,  efficient  and  versatile.  In  other  words,  it  must 
consistently  provide  accurate  results  for  a  reasonable  cost  for  a  wide  range  of  problems.  The  first 
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Figure  8.  Convergence  of  Residual  Norm  and  Reflection  Coefficient 
using  Conjugate  Gradient  Method 

criterion,  effectiveness,  will  be  established  by  the  validation  cases  to  be  discussed  in  Section  6.  It  is  the 
result  of  a  faithful  implementation  of  the  formulas  from  the  preceding  sections.  The  second  criterion, 
efficiency,  is  partly  a  consequence  of  inherent  properties  of  the  solution  method,  but  also  results  from 
design  of  algorithms  to  minimize  storage  and  computation  time.  Those  algorithms  are  discussed  in 
Sections  5.1  and  5.2.  The  last  criterion,  versatility,  is  also  partly  a  consequence  of  the  method,  but  it 
cannot  be  realized  without  a  means  of  defining  complex  geometries.  That  problem  was  solved  by 
interfacing  the  code  to  commercial  computer-aided  design  software,  as  discussed  in  Section  5.3. 

5.1.  Computational  Algorithms 

The  overall  design  for  the  computer  program  TWOPORT  is  shown  in  Figure  9,  most  of  which 
is  self-explanatory.  The  main  loop  may  be  executed  repeatedly  for  increments  in  frequency.  The  major 
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Figure  9.  Flow  Diagram  for  Program  TWOPORT 

subroutines  involved  in  each  step  are  annotated  on  the  right.  (The  program  is  written  in  ANSI  standard 
FORTRAN  and  has  been  executed  on  a  VAX®  8650  minicomputer  and  a  Cray  Y-MP®.) 

The  serious  work  begins  in  the  second  block,  which  imports  a  geometry  description  that  was 
produced  by  the  CAD  software.  That  description,  illustrated  in  Figure  10,  consists  of  a  tetrahedron  mesh 
defined  by:  (a)  a  list  of  vertices  with  their  x,y.z  coordinates,  and  flags  that  indicate  whether  they  are  on 
an  aperture  and/or  on  one  or  more  conducting  surfaces  (allowing  for  the  possibility  that  a  node  may  be 
at  the  intersection  of  conducting  surfaces);  (b)  a  list  of  tetrahedra  with  the  four  nodes  defining  each  and 
the  index  of  the  material  filling  it;  and  (c)  a  list  of  the  materials  with  their  permittivity  and  permeability. 
The  geometry  conversion  is  an  additional  step  that  is  required  to  accommodate  the  edge-based  vector 
expansion  functions:  (a)  create  a  list  of  edges  with  their  terminal  node  indices;  and  (b)  create  a  new  list 
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Figure  10.  Example  Input  Geometry  File  (as  Created  by  Translator 
Program  from  CAD  Universal  File;  italics  are  annotation) 

of  tetrahedra  with  the  six  edges  that  define  each  one.  Note  that  during  this  process,  if  two  connected 
nodes  are  on  the  same  conductor,  then  that  edge  is  not  included  in  the  edge  list,  so  it  will  be  ignored  in 
subsequent  calculations.  This  mechanism  forces  the  field  along  those  edges  to  be  zero,  thereby  enforcing 
conductor  boundary  conditions.  However,  when  two  nodes  are  on  different  conductors,  there  may  be 
a  field  component  along  the  edge  linking  them,  and  so  that  edge  must  be  accounted  for.  That  is  the 
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reason  why  the  geometry  model  must  identify  which  specific  conductors  a  node  is  on,  an  not  simply 
whether  or  not  it  is  on  one. 

The  interior  matrix  computations  are  Eqs.(34)  and  (36).  These  are  performed  one  cell  at  a  time, 
as  depicted  in  Figure  1 1 . 

When  using  row-column  storage  in  the  LUD  version  of  the  program,  the  waveguide  contributions 
are  not  calculated  as  separate  matrices  as  implied  by  Figure  9.  They  are  merely  added  to  the  correct 
locations  of  the  matrix  that  already  contains  the  interior  (FEM)  results.  This  conserves  storage  by 
requiring  only  one  matrix  instead  of  three.  As  Figure  12  shows,  these  calculations  are  carried  out  one 
mode  at  a  time  for  m/i=  10,1 1,...  For  each  mode  the  inner  products  are  generated  for  each  edge 
and  each  polarization  p.q.  Then,  in  the  last  block,  the  matrix  elements  are  updated  for  each  pair  r,s  in 
the  aperture.  An  important  detail  that  is  not  reflected  in  Figure  10  is  that  each  (the  dominant  mode) 
is  saved  for  later  use  in  computing  the  reflection  and  transmission  coefficients.  A  single  array  is  used 
for  both  since  edges  in  ports  1  and  2  are  mutually  exclusive,  so  their  nonzero  entries  do  not  overlap. 

5.2.  Sparse  Storage  Structure  and  Conjugate  Gradient  Algorithm 

The  first  version  of  program  TWOPORT  used  an  ordinary  row-column  storage  format.  Relative¬ 
ly  small-scale  problems  can  thus  be  solved  using  LU  decomposition  by  standard  library  routines.  It  also 
allowed  us  to  gain  confidence  in  the  code  without  the  added  uncertainty  of  matrix  solver  convergence. 
The  LUD  version  is  still  recommended  for  problems  having  1000  edges  or  less. 

An  existing  CGM  solver  |14)  was  used  in  conjunction  with  the  row-column  storage  to  gather 
convergence  data  for  some  small  test  problems.  However,  taking  full  advantage  of  the  CGM  requires 
using  a  sparse  storage  format  tailored  to  the  form  of  matrix  generated  by  this  code.  The  structure  that 
was  adopted  uses  three  separate  arrays:  one  for  the  interior  (sparse)  matrix  generated  by  the  finite 
element  computations;  and  one  each  for  the  waveguide  apertures,  which  are  dense.  The  dense  matrices 
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Figure  1 1 .  Structure  of  Finite  Element  Interior  Computations 
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Figure  12.  Structure  for  Waveguide  Computations 


29 


are  stored  in  row-column  form. 


The  interior  matrix  entries  are  stored  in  a  column  array  in  the  order  in  which  they  are  created. 
That  order  depends  only  on  the  sequence  in  which  edge  pairs  are  encountered.  For  instance,  if  cell  #1 
includes  edges  3,4,  and  7  (although  a  tetrahedron  has  six  edges,  those  that  are  on  perfect  conductors  are 
not  included  in  the  edge  list,  and  so  a  cell  may  have  fewer  than  6  numbered  edges),  then  the  first  entries 
to  be  created  are  ^33,  534,  S37,  ^43,  ^44.  547,  ^73,  574,  S77,  and  these  are  stored  in  A(l)  -  A(9).  If  any 
of  these  pairs  are  shared  by  other  cells,  then  those  cells’  contributions  are  added  to  the  existing  entries 
in  A.  Two  integer  arrays  lA  and  JA  keep  track  of  the  row  and  column  locations,  respectively,  of  the 
entries  in  A,  that  is;  1A  =  (3,3,3.4,4,4,7,7,7,...|  and  JA  =  |3,4,7,3,4,7,3,4.7....1.  To  aid  in  searching  for 
existing  entries,  two  arrays  are  maintained  to  point  to  the  first  and  last  entries  for  each  row.  for  example; 
NF(3)  =1,  NL(3)  =  3,  NF(4)  =  4,  NL(4)=6.  etc.  Since  calculations  for  subsequent  cells  will  involve 
searching  A  to  determine  whether  a  given  edge  pair  already  exists,  these  arrays  NF  and  NL  speed  the 
process  by  bounding  the  search. 

Implementing  Eqs.(44)-(49)  requires  approximately  2MN  multiplies,  where  N  is  the  number  of 
edges  and  M  is  the  average  edge  connectivity  {MN  is  the  total  number  of  nonzero  entries  in  A).  M  is 
independent  of  the  problem  size,  and  is  usually  around  10.  If  convergence  is  obtained  in  N  iterations, 
then  the  execution  time  is  OQMN-)  operations.  On  the  other  hand,  LU  decomposition  requires  0(A  ) 
operations. 

The  algorithm  as  programmed  is  essentially  a  straightforward  implementation  of  Eqs.(43)-(49). 
The  initial  guess  for  the  solution  is  the  incident  field  for  the  first  frequency.  For  subsequent  frequencies, 
the  converged  result  from  the  previous  frequency  is  used  as  the  initial  Aq. 

5.3.  I-DEAS^"^  Interface 

The  commercial  software  I-DEAS^^  is  a  system  of  programs  oriented  toward  structural  analysis 
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[15].  Only  a  few  of  its  components  are  used  in  this  work,  most  notably,  the  solid  object  modeling  and 
mesh  generation.  The  "solid  object"  is  the  region  denoted  B.  It  is  the  entire  ccvity,  excluding  the 
interiors  of  conducting  obstacles.  For  one  test  case  the  cavity  was  simply  a  short  section  of  coaxial  line, 
for  which  the  solid  object  is  the  tube  shown  in  Figure  13a. 

Mesh  generation  is  an  automated  process  that  creates  nodes  on  the  boundaries  and  in  the  interior 
of  the  solid  object  so  that  the  distance  between  any  node  to  its  nearest  neighbor  is  less  than  or  equal  to 
a  user-specified  global  element  edge  length.  Then  the  nodes  are  grouped  into  tetrahedral  elements 
bounded  by  groups  of  four  nodes.  The  result  is  shown  in  Figure  13b. 

There  are  a  few  preliminary  steps  before  the  automatic  mesh  generation  may  be  executed.  First, 
points  and  curves  are  created  from  the  solid  object  in  the  "geometry  definition"  task.  Then  the  curves 
are  linked  together  to  form  mesh  areas  in  the  "mesh  generation"  task.  Last,  the  mesh  areas  are  assembled 
into  mesh  volumes. 


(b) 

(a)  Solid  Object;  (b)  Tetrahedron  Mesh 
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The  electromagnetic  analysis  code  will  need  to  have  nodes  identified  by:  (a)  the  port  number,  if 
they  are  in  an  aperture;  and  (b)  conducting  surface  numbers(s)  if  applicable.  The  port  identifier  is 
conveyed  through  a  boundary  condition  variable;  the  z  component  of  "nodal  force,"  set  to  1  or  2  for 
waveguide  A  or  B,  respectively.  The  translator  program  (TRANSL8)  will  examine  that  variable  in  the 
I-DEAS  universal  file  to  determine  which  nodes  are  waveguide  aperture  nodes.  The  conducting  surfaces 
are  defined  by  grouping  nodes  into  "permanent  groups."  The  translator  program  looks  for  permanent 
groups  in  the  universal  file  and  uses  those  to  identify  each  node’s  inclusion  in  up  to  six  conductor 
surfaces.  Last,  the  material  properties  in  the  cells  must  be  identified.  I-DEAS,  being  oriented  to  me¬ 
chanics,  does  not  provide  menu  options  for  defining  permittivity  and  permeability.  However,  the  "physi¬ 
cal  properties"  table  has  nine  general-purpose  data  fields  called  a  "keyopt  array,"  and  the  first  four  of 
these  are  used  for  Re(e^),  Im(£,),  ReOi^)  and  Im(/if).  Any  number  of  dielectric  regions  may  be  defined 
by  grouping  elements  and  assigning  them  to  different  physical  properties  tables. 

The  last  step  inside  I-DEAS  is  to  write  the  universal  file.  That  file  contains  a  great  deal  of 
information  not  needed  by  the  electromagnetics  code.  The  translator  program  filters  out  the  superfluous 
information  and  writes  an  output  file  in  a  form  (see  Figure  10)  that  is  easier  to  inspect  and  modify. 

6.  VALIDATION  AND  DEMONSTRATION 

This  section  presents  test  results  for  the  computer  program  TWOPORT.  The  tests  were  intended 
to  serve  three  purposes:  first,  to  verify  that  the  code  produces  correct  results;  second,  to  show  that  it  is 
useful  for  problems  that  have  practical  applications;  and  third,  to  assess  the  modeling  requirements  such 
as  mesh  density  and  number  of  waveguide  modes. 
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6.1.  Waveguides  and  Waveguide  Discontinuities 


6.1.1.  WAVEGUIDE  PROPh»^  .- ION. 

Referring  again  to  Figure  1 ,  suppose  that  the  inlet  and  outlet  waveguides  are  identical  in  size  and 
shape,  and  that  the  "cavity"  is  another  section  of  the  same  waveguide  with  length  d.  Then  the  magnitudes 
of  Sji  and  $21  should  be  0.0  and  1.0,  respectively,  since  there  are  discontinuities  or  lossy  materials. 
Furthermore,  the  phase  of  Soj  should  be  260°dl\^,  where  Xg  is  the  guide  wavelength  of  the  dominant 
mode. 

The  coaxial  waveguide  section  shown  earlier  in  Figure  13  has  the  dimensions  of  an  APC-7mm 
connector  (a=  1 .502  mm,  b  =  3 .5  nun)  and  it  is  4  mm  long.  The  guide  wavelength  is  the  same  as  the  free 
rpace  wavelength  since  there  is  no  dielectric  and  the  dominant  mode  is  TEM.  Calculations  by  TWO- 
PORT  over  a  range  of  frequencies  are  shown  in  Figure  14.  The  agreement  with  the  theory  is  quite  good. 


Figure  14.  Propagation  Phase  Through  4  mm  APC-7  Coaxial  Waveguide  Section 
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indicating  that  the  finite  element  mesh  does  not  introduce  any  dispersion.  At  all  frequencies  considered, 
ISji  1  and  IS21I  were  within  lx  10^  of  the  correct  values. 

Figure  15  is  the  mesh  for  a  section  of  standard  X-band  rectangular  waveguide  (3=22.86  mm, 
b=  10.16  mm)  with  length  d  =  b/2  =  5.08  mm.  The  S21  phase  computed  by  TWOPORT,  shown  in  Figure 
16.  agrees  very  well  with  the  expected  results  for  the  TEjo  mode  given  by 


=  -360“ 


(50) 


For  each  frequency  the  magnitude  of  S,,  was  less  than  0.01  and  that  of  S21  was  greater  than  0.9999. 
The  rectangular  waveguide  results  are  important  because  they  demonstrate  that  the  finite  element  imple¬ 
mentation  is  correctly  enforcing  conductor  boundary  conditions  at  interior  corners. 

Figure  17  shows  the  mesh  for  a  section  of  circular  waveguide.  The  mesh  generator  placed  the 


Figure  15.  Finite  Element  Mesh  for  a  Rectangular  Waveguide  Section 
(a=22.86  mm,  b=  10.16  mm,  d=5.08  mm) 
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Figure  16.  Propagation  Phase  Delay,  Rectangular  Waveguide  Section 


Figure  17.  Finite  Element  Mesh  for  a  Circular  Waveguide  Section 
(radius =0.5  m,  length =0.4  m) 
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Figure  18.  Propagation  Phase  Through  Circular  Waveguide  Section 
(single  mode  propagation  range  is  176-231  MHz) 

nodes  on  the  surface  of  the  cylinder,  and  consequently  the  collection  of  tetrahedra  contains  slightly  less 
total  volume  than  the  original  cylinder.  It  is  necessary  to  scale  the  mesh  in  x  and  y  to  correct  for  that 
modeling  error.  (A  scaling  option  is  incorporated  in  the  translator  program.)  As  Figure  18  shows,  the 
calculated  and  predicted  phase  of  the  dominant  TEj  j  mode  agree  quite  well  when  this  volume  correction 
is  used. 

6.1.2.  WAVEGUIDE  IRISES. 

Thin  metallic  obstructions  oriented  transverse  to  the  waveguide  axis  have  applications  as  elements 
of  microwave  circuits,  especially  filters  [16].  Our  interest  in  them  lies  in  their  tendency  to  generate 
higher-order  waveguide  modes,  making  them  useful  as  test  cases  for  validating  the  aperture  interactions 
through  those  higher  order  modes. 

The  first  two  iris  test  cases  used  a  short  section  of  APC-7  coaxial  waveguide,  just  as  in  Figure 
13,  except  that  now  part  of  the  inlet  aperture  will  be  blocked.  Network  analyzer  measurements  were 
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made  with  foil  obstructions  placed  between  coaxial  adapters,  blocking  either  one  quarter  (a  quadrant)  or 
one  half  of  their  common  aperture.  Figures  19a  and  19b  show  comparisons  between  measured  and 
calculated  S||,  both  magnitude  and  phase,  for  these  test  cases. 

The  third  iris  test  case  was  a  diagonal  obstruction  in  X-band  rectangular  waveguide.  The  finite 
element  mesh  model  is  shown  in  Figure  20a.  Figure  20b  shows  the  mesh  with  cells  above  z=5  mm 
removed,  and  shading  added  to  identify  the  iris  shape  and  location.  Figure  21  is  a  comparison  of  the 
TWOPORT  calculations  and  network  analyzer  measurements  of  Su. 

The  close  agreement  between  the  measurements  and  calculations  for  these  waveguide  iris  cases 
serves  to  validate  the  higher-order  mode  computations,  which  are  fairly  extensive,  especially  for  the 
coaxial  waveguide.  They  include  calculations  of:  (a)  cutoff  wavelength;  (b)  mode  functions;  and  (c) 
integration  of  inner  products  of  mode  functions  and  finite  elements.  The  circular  waveguide  mode 
functions  were  validated  by  a  step  discontinuity  problem,  discussed  in  the  following  section. 

6.1.3.  CIRCULAR  WAVEGUIDE  MODE  CONVERTER. 

A  step  discontinuity  in  a  circular  waveguide  is  a  fairly  simple  mode  converter  (17).  The  outlet 
waveguide  is  slightly  larger  than  the  inlet  waveguide  so  that  it  supports  the  TM|]  mode  as  well  as  the 
TE,,  mode.  It  is  the  abrupt  change  in  radius  that  excites  the  TMj,  mode. 

Masterman  and  Clarricoats  1 17]  give  comparisons  of  multimode  calculations  against  measured  data 
for  test  cases  in  which  the  outlet  waveguide  has  radius  of  35.56  mm,  and  several  different  inlet  waveguide 
radii.  Their  results  are  presented  in  terms  of  a  "mode  conversion  ratio,"  M\ 

M  .  (sr) 

1  D^h^^(J3=a)\ 

D]  and  D2  are  the  excitation  coefficients  for  the  TEu  and  TM,]  modes,  respectively,  in  the  outlet  wave¬ 
guide,  and  and  hp2  are  the  radial  components  of  their  orthonormal  mode  functions.  Hence  M  is  the 
relative  strength  of  the  two  modes  measured  at  the  waveguide  wall,  p  =  a.  The  ratio  of  hp2  to  is 
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Figure  19.  Measured  and  Calculated  Reflection  Coefficient  vs.  Frequency 
for  Coaxial  Irises  (APC-7mni  coax):  (a)  180°  Opening;  (b)  270°  Opening 


REFLECTION  COEFFICIENT  PHASE  (deg.)  ^  REFLECTION  COEFFICIENT  PHASE  (deg.) 


Figure  20.  Finite  Element  Mesh  for  Diagonal  Iris  in  Rectangular  Wave¬ 
guide;  (a)  Complete  Mesh;  (b)  Cells  above  z—5  mm  removed 


Figure  21.  Measured  and  Calculated  Sj|  for  Diagonal  Iris  in  Rectangular 
Waveguide  (X-band  waveguide) 


independent  of  frequency,  and  is  3.79  dB  for  the  35.56  mm  radius  waveguide. 

Figure  22  shows  the  finite  element  meshes  used  for  two  test  cases,  each  with  a  different  inlet 
waveguide  radius;  (a)  26.67  mm  (1.05  in.);  and  (b)  29.21  mm  (1.15  in.).  The  meshes  produced  by  the 
CAD  software  were  scaled  appropriately  to  ensure  that  the  included  volume  was  correct.  In  both  models, 
the  total  height  in  z  is  10  mm,  evenly  divided  between  the  two  waveguide  sections. 

The  excitation  coefficients  for  the  higher  order  modes  are  computed  by  the  TWOPORT  code,  and 
from  these  it  is  a  simple  matter  to  find  the  mode  conversion  ratio.  The  results,  shown  in  Figure  23, 
agree  with  the  multimode  calculations  to  within  1  dB  over  the  entire  frequency  range.  These  results 
demonstrate  a  flexibility  of  our  approach  and  the  TWOPORT  code  for  analyzing  waveguide  structures 
that  may  have  many  propagating  modes. 

6.2.  Waveguide  Transitions 

Transitions  between  different  types  of  waveguiding  media  are  an  important  class  of  devices  that 
can  be  modeled  using  the  hybrid  finite  element  method.  The  circular  waveguide  mode  converter  was  a 
simple  case.  This  section  considers  two  structures  that  are  somewhat  more  complex:  a  coax -stripline 
transition;  and  a  coaxial  end-wall  launcher  for  rectangular  waveguide. 

6.2.1.  COAX-STRIPLINE  TRANSITION. 

Figure  24a  shows  a  short  section  of  stripline  (with  air  dielectric)  enclosed  in  a  rectangular 
conducting  box.  The  openings  at  each  end  are  the  appropriate  size  for  SMA  coaxial  adapters.  The  mesh 
geometry  in  Figure  24b  shows  only  the  cells  below  the  stripline.  Shading  has  been  added  to  identify 
those  cells  that  border  on  the  strip  conductor,  the  coax  dielectric,  and  end  of  the  coax  center  conductor, 
which  is  trimmed  flat  where  it  meets  the  stripline.  The  cross-sectional  dimensions  for  both  the  coax  and 
the  stripline  are  given  in  Figure  24c.  The  stripline  width  was  chosen  for  a  characteristic  impedance  of 
50  n,  using  formulas  from  Gupta  et.  al.  (Reference  18,  pp.  57-60).  The  computed  reflection  and  trans- 
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MODE  CONVERSION  RATIO  (dB) 


Figure  22.  Finite  Element  Meshes  for  Circular  Waveguide  Step  Discontinuities 
(Mode  Converter),  with  32=  1 .4  in.;  (a)  a]  =  1 .05  in.;  (b)  a,  =  1 . 15  in. 
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Figure  23.  Mode  Conversion  Ratio  for  Circular  Waveguide  Mode  Converters, 
Hybrid  Finite  Element  (HFEM)  and  Multimode  |17]  Calculations 
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Figure  24.  Stripline  Test  Case:  (a)  Complete  Problem  Geometry;  (b)  Lower 
Half  of  Finite  Element  Mesh;  and  (c)  Dimensions  at  Coax/Stripline  Junction 


mission  coefficients  for  a  4  mm  long  stripiine  are  shown  in  Figure  25.  As  expected,  the  reflection 
coefficient  is  small  at  the  lower  frequencies  because  the  coax  and  the  stripline  have  the  same  characteristic 
impedance.  But  the  junction’s  reactances  become  more  significant  with  increasing  frequency,  as  indicated 
by  a  gradual  rise  in  |S|,  j,  until  it  decreases  again  approaching  18.75  GHz,  where  there  is  a  resonant 
cancellation  because  the  stripiine  is  one  quarter  wavelength  long  at  that  frequency.  These  results  are 
important  because  they  demonstrate  that  conductor  boundary  conditions  are  correctly  enforced  at  sharp 
conductor  edges  and  at  multi-conductor  Junctions. 

6.2.2.  COAX-RECTANGULAR  WAVEGUIDE  LAUNCHER. 

A  design  described  by  Tang  &  Wong  1 19|  for  an  end-wall  transition  between  coax  and  rectan¬ 
gular  waveguide  is  shown  in  Figure  26a.  Their  design  used  X-band  rectangular  waveguide  (10  mm  x 
23  mm).  The  coax  inner  and  outer  radii  were  1  mm  and  2.3  mm,  respectively.  The  coax  port  was 
centered  in  the  y-direction  (the  waveguide  height)  but  it  was  slightly  off  center  in  the  x  direction,  by 
1/lOth  the  waveguide  width.  The  2  mm-thick  post  was  centered  8.4  mm  from  the  end  wall.  The  probe 
(an  extension  of  the  coax  center  conductor)  was  9.4  mm  long  overall.  For  convenience  of  modeling,  the 
round  conductors  are  modeled  as  eight-sided  polygons  and  the  corner  between  the  probe  and  post  was 
sliced  at  a  45°. 

Figure  26b  is  a  cutaway  view  of  the  finite  element  mesh  for  this  object.  It  was  fairly  challenging 
in  terms  of  the  geometry  generation  requirements,  and  the  success  is  attributable  to  the  flexibility  of  the 
CAD  software.  The  Figure  illustrates  that  the  interiors  of  solid  conductors  are  represented  as  voids  in 
the  finite  element  mesh.  As  usual,  the  nodes  that  lie  on  those  surfaces  were  appropriately  tagged  to 
identify  them  for  the  electromagnetic  analysis  code. 

Figure  27  is  the  calculated  reflection  coefficient  vs.  frequency  for  this  launcher  design.  As 
expected,  it  performs  well  (less  than  2:1  VSWR)  over  most  of  the  X-band  (8  -  12  GHz). 
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REFLECTION  COEFFICIENT  MAGNITUD 


6.3.  Printed  Circuit  Devices 


Up  to  now  this  report  has  only  discussed  structures  with  air  dielectric.  Yet  one  of  the  most 
important  capabilities  of  the  finite  element  method  is  its  ability  to  deal  with  inhomogeneous  dielectrics. 
This  section  will  discuss  results  for  two  microstrip  structures:  first,  a  calculation  of  guide  wavelength 
in  a  microstrip  line;  and  second,  calculations  of  S  parameters  for  a  microstrip  meander  line. 

6.3.1.  MICROSTRIP  TRANSMISSION  LINE. 

Figure  28  shows  the  mesh  for  a  length  of  75  fim  wide  microstrip  line  on  a  100  /xm-thick  substrate 
(Gallium  Arsenide,  €^=12.9).  The  coaxial  connector  dimensions  are  a  =  43  /xm  and  b  =  100  ^m, 
chosen  for  convenience  and  an  impedance  match  to  the  50  Q  microstrip.  The  substrate  is  1000  /xm  wide 
and  the  cavity  height  is  500  ^m  total,  with  400  ^m  of  air  between  the  microstrip  and  the  cavity’s  top 
wall.  These  dimensions  ensure  that  there  is  a  space  of  4-5  line  widths  between  the  microstrip  and  the 
enclosure,  which  is  adequate  to  ensure  that  the  package  does  not  influence  the  guide  wavelength  or 
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CASE  OUTUNE 


Figure  28.  Finite  Element  Mesh  for  Microstrip  Line  with  Coaxial  Connectors 
characteristic  impedance  [20]. 

Determining  the  properties  of  the  microstrip  by  itself  requires  calibrating  out,  or  "de-imbedding" 
the  effects  of  the  coaxial  connectors.  Using  a  scaling  option  in  the  TWOPORT  code,  the  S  parameters 
were  computed  first  for  a  line  length  of  500  /xra;  then  second  with  the  same  finite  element  mesh  stretched 
in  z  to  a  length  of  750  /xm.  The  differential  phase  in  $21  was  35°  for  the  test  frequency  of  40  GHz.  This 
yields  2.57  mm  for  the  guide  wavelength,  while  quasi-static  formulas  developed  by  Wheeler  [21]  give 
2.54  mm.  This  close  agreement  indicates  that  the  model  is  correctly  predicting  the  fields  at  interfaces 
between  two  dielectrics,  even  in  the  presence  of  sharp  conducting  edges. 

6.3.2.  MICROSTRIP  MEANDER  LINE. 

Since  the  modeling  method  accurately  accounts  for  all  relevant  boundary  conditions  for  a  simple 
microstrip  line,  it  should  be  obvious  that  it  can  correctly  predict  the  performance  of  most  passive  printed- 
circuit  RF  components.  A  microstrip  meander  line  was  chosen  as  a  demonstration  case  since  measure¬ 
ments  were  available.  The  circuit  dimensions  (in  /im)  and  the  finite  element  mesh  (substrate  cells  only) 
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are  shown  in  Figure  29. 

The  measurement  reference  planes  are  at  the  centers  of  the  pads  at  each  end  of  the  circuit,  while 
those  used  in  the  calculations  are  at  the  points  where  the  lines  begin  to  taper  from  75  ^m  down  to  25  ^m, 
so  the  difference  must  be  accounted  for  when  comparing  the  d?ta.  The  circuit  layout  shows  the  location 
of  "via  holes"  that  proviue  a  ground  for  the  measurement  probes.  At  the  contact  points,  the  75  nm  center 
conductor  and  50  /xm-wide  slots  form  a  50  0  coplanar  waveguide.  Using  formulas  given  by  Rowe  and 
Lao  (22j,  the  effective  relative  permittivity  for  that  transmission  line  structure  was  found  as  €gff=6.29. 
The  measured  Soj  phase  was  corrected  by  adding  360°  *(75  The  calculated  S21  phase 

was  also  corrected  by  subtracting  the  excess  phase  introduced  by  the  coaxial  connectors,  computed  from 
a  straight  length  of  microstrip  line  as  discussed  in  the  preceding  section.  Figure  30  shows  that  the 
measurements  and  calculations  agree  very  well  over  the  frequency  range  from  1  to  26  GHz.  Note  that 
the  measurement  apparatus  uses  coplanar  waveguide  probes,  while  the  calculation  simulates  coaxial 


Figure  29.  Microstrip  Meander  Line;  (a)  Wafer  Metallization  Dimensions 
(in  fim);  and  (b)  Finite  Element  Mesh  for  Substrate 
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Figure  30.  Comparison  of  Measured  and  Calculated  S21  Phase  for 
Microstrip  Meander  Line 

connectors.  The  CPW-microstrip  and  coax-microstrip  transitions  have  different  reactances,  which 
accounts  for  the  difference  in  the  phase  slope  between  measurement  and  calculation  in  Figure  30. 

This  concludes  discussion  of  the  validation  results.  Together  they  demonstrate  the  effectiveness, 
reliability  and  versatility  of  the  hybrid  finite  element  method.  The  final  part  of  this  section  presents 
conclusions  that  may  be  drawn  from  these  tests  regarding  rule-of-thumb  modeling  requirements  such  as 
number  of  waveguide  modes,  number  of  finite  elements  and  typical  execution  time. 

6.4.  Modeling  Requirements 

6.4.1.  HIGHER-ORDER  MODES. 

All  of  the  validation  cases  used  many  more  waveguide  modes  than  were  actually  required.  The 
additional  matrix  fill  calculations  are  not  very  significant  for  rectangular  waveguides,  but  they  can  be 
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somewhat  expensive  for  coax.  If  computation  time  is  critical,  it  may  be  important  to  know  how  many 
modes  are  actually  needed  to  obtain  an  accurate  answer. 

Approximate  matrix  fill  execution  times  are  as  follows  for  the  VAX  8650  computer,  given  in 
terms  of  CPU  seconds  per  mode  per  aperture  edge:  rectangular,  2.5  ms;  circular,  25  ms;  and  coaxial, 
250  ms.  The  circular  and  coaxial  mode  function  integrations  both  use  a  16-point  quadrature  rule, 
meaning  that  there  are  16  integration  sample  points  in  each  aperture  triangle. 

The  TWOPORT  code  includes  a  calculation  of  the  excitation  coefficients  for  all  of  the  modes 
included  in  the  matrix  calculation.  These  coefficients  are  the  C,’s  from  Eq.(7).  Figure  31  shows  the 
results  for  the  two  coaxial  irises  discussed  in  Section  6.1.2.  Apparently,  there  is  a  "noise  floor"  at  -30 
to  -40  dB  due  to  gridding  (discretization)  effects.  These  iris  discontinuities  evidently  only  excite  TE 
modes,  and  only  higher  orders  of  m.  This  is  consistent  with  expectations,  since  higher  orders  of  n 
correspond  to  functions  with  more  cycles  in  the  radial  direction  and  should  not  be  excited  by  radially- 
invariant  obstructions.  Since  coaxial  and  circular  modes  degenerate  into  sinfmif))  and  cos(m<f>)  components 
(see  Appendix  B),  both  must  be  included  to  represent  the  fields  at  arbitrary  obstructions,  which  may  not 
be  aligned  with  the  coordinate  axes. 

Evidently,  only  modes  /n/i=01,ll,  and  21  are  needed  for  these  coaxial  irises.  This  conclusion 
is  reinforced  by  Figure  32,  which  shows  the  computed  values  of  |S]|  |  for  different  numbers  of  modes. 
Adding  modes  higher  than  21  did  not  produce  any  change  in  the  reflection  coefficient,  so  they  are  not 
shown  in  Figure  32. 

The  coax -stripline  transition  discussed  in  Section  6.2.1  is,  by  contrast,  a  longitudinal  (parallel  to 
z)  rather  than  a  lateral  (transverse  to  z)  discontinuity.  Consequently,  it  excites  higher  order  TM  modes. 
Its  excitation  coefficients  are  shown  in  Figure  33.  The  strength  of  these  modes  is  abnormally  high 
because  the  coax  diameter  is  so  much  larger  than  the  cavity  height,  so  these  results  are  not  representative 
of  coaxial  transitions  to  stripline  or  microstrip.  In  fact,  the  microstrip  circuits  in  Sections  6.3.1  and  6.3.2 
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Figure  32.  Reflection  Coefficient  Magnitude  vs.  Frequency  Calculated 
with  Varying  Numbers  of  Higher  Order  Modes,  270°  Coaxial  Iris 
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Figure  33.  Higher-Order  Coaxial  Mode  Excitation  Coefficients  at  Coax- 

Stripline  Transition 


did  not  generate  any  appreciable  higher  order  modes  because  the  coax  outer  radius  was  the  same  as  the 
substrate  height. 

The  last  issue  to  be  addressed  is  whether  there  is  in  fact  any  benefit  to  using  higher -order  modes. 
The  alternative  is,  of  course,  to  extend  the  finite  element  mesh  into  the  connecting  waveguides  far  enough 
that  any  higher  modes  excited  by  the  cavity’s  contents  and/or  apertures  have  decayed  to  insignificance 
[23).  The  answer  depends  on  two  factors:  first,  the  mode  excitation  strengths;  and  second,  their  atten¬ 
uation  constants  in  the  waveguides,  which  depend  mainly  on  the  ratio  of  the  frequency  to  the  mode  cutoff 
frequency.  Consider  the  following  two  extreme  examples;  First,  the  microstrip  circuits  do  not  generate 
higher-order  modes  in  the  coaxial  waveguides.  Therefore,  the  mesh  may  be  terminated  at  the  coax 
apertures  even  when  only  the  dominant  mode  is  used.  At  the  other  extreme,  the  one-quadrant  coaxial 
iris  (Figure  30b)  excites  the  TE,i  mode  at  about  -10  dB.  Its  cutoff  frequency  in  APC-7  coax  is  about 
20  GHz.  At  18  GHz  it  decays  by  29  dB  per  wavelength,  and  for  it  to  be  negligible  (below  -30  dB),  the 
mesh  would  need  to  extend  roughly  %X  into  the  waveguide  in  each  direction.  The  finite  element  mesh 
and  the  interaction  matrix  would  then  contain  an  additional  500-1500  edges  and  unknowns,  respectively. 
Thus  it  is  clear  that  in  some  cases  at  least,  the  use  of  higher-order  modes  is  very  beneficial. 

6.4.2.  DISCRETIZATION  REQUIREMENTS. 

The  discretization  for  most  of  th  ♦est  cases  is  summarized  in  Table  1.  It  lists  tfie  total  volume 
in  cubic  wavelengths  at  the  highest  frequency,  the  total  number  of  cells  and  the  number  of  cells  per  cubic 
wavelength.  Accurate  results  for  S  parameters  may  be  obtained  with  as  few  as  10  nodes  or  edges  per 
wavelength  in  each  linear  direction.  Finer  sampling  may  be  required  in  some  instances.  For  example, 
the  circular  waveguide  mesh  in  Figure  13  was  not  adequate  for  the  mode  converter  problem.  It  produced 
fairly  good  answers  for  the  case  with  the  smaller  inlet  waveguide  radius,  but  its  results  were  off  by  as 
much  as  4  dB  for  the  other  case.  The  reason  is  that  the  number  of  integration  points  over  the  collection 
of  triangles  in  the  aperture  was  not  enough  to  give  an  accurate  result  for  inner  products  with  the  higher- 
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order  mode  functions.  No  more  than  20  nodes  per  wavelength  in  any  linear  direction  should  be  required, 
but  in  dense  dielectrics,  such  as  GaAs,  the  wavelength  is  much  smaller  than  in  free  space.  In  many 
instances  the  mesh  density  was  dominated  by  fine  details  in  the  geometry.  I-DEAS  attempts  to  make  each 
tetrahedron  as  nearly  equilateral  as  possible,  which  does  not  allow  the  node  spacing  to  "relax"  very 
abruptly.  The  rationale  is  that  highly  elongated  cells  may  cause  ill  conditioning  in  certain  finite  element 
applications.  The  consequence  here  is  that  there  are  many  more  cells  than  are  required  to  make  accurate 
field  computations. 

6.4.3.  MATRIX  SOLUTION. 

Figure  34  summarizes  the  execution  times  observed  for  the  various  test  cases.  All  of  these  were 
on  the  VAX  8650  using  single  precision.  The  LU  decomposition  was  performed  by  the  IMSL  routine 
L2LCG,  which  is  optimized  to  use  a  vector  floating  point  coprocessor.  LU  decomposition  was  also 
executed  on  a  Cray  Y-MP,  and  typically  ran  about  60  times  as  fast  as  it  did  on  the  VAX. 


100  1000  10000  100000 
NUMBER  OF  UNKNOWNS 


Figure  34.  Execution  Time  vs  Number  of  Unknowns  for  LU  Decomposition 
and  Conjugate  Gradient  Method  (VAX  8650) 
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Table  1.  Summary  of  Finite  Element  Discretization  for  Test  Cases 


Test  Case 

Reference 
Figure  No. 

Volume 

No.  Cells 
in  Mesh 

Cells 

ARC -7  Coax 

13 

0.0135 

847 

63x10^ 

Rectangular 

Waveguide 

15 

0.149 

308 

2067 

Circular 

Waveguide 

17 

0.120 

342 

2850 

Circular  W/G 

Mode  Converter 

22a 

0.376 

1504 

4000 

Stripline 

24b 

0.0138 

506 

37x10^ 

Launcher 

26 

0.2944 

2442 

8294 

Microstrip  Line 

27 

5.7x10^ 

4662 

8.2x10^ 

Meander  Line 

28 

9.4x10^ 

4190 

4.5x10'^ 

The  conjugate  gradient  solver  is  not  vectorized,  yet  for  small  numbers  of  unknowns  it  runs  almost 
as  fast  as  LU  decomposition.  It  is  clearly  a  better  choice  for  large  problems,  but  the  scattering  of  points 
in  Figure  34  reflects  the  disadvantage  that  the  solution  time  is  difficult  to  predict  a  priori. 

1.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  hybrid  finite  element/waveguide  mode  approach  is  clearly  an  effective  tool  for  modeling  the 
electromagnetic  fields  in  a  wide  variety  of  radio  frequency  devices.  The  success  of  the  implementation 
presented  in  this  report  indicates  that  a  number  of  features  of  this  approach  are  appropriate  and  correct, 
specifically:  (a)  the  representation  of  the  boundary  value  problem  as  a  weak  form  of  the  Helmoltz 
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equation  for  the  electric  field;  (b)  the  enforcement  of  transverse  magnetic  field  continuity  across  wave¬ 
guide  apertures  using  the  boundary  term  of  the  functional;  (c)  the  representation  of  waveguide  fields  by 
mode  sums  and  the  implementation  of  the  higher-order  mode  functions;  (d)  the  enforcement  of  the 
divergence  condition  and  conductor  boundary  conditions  using  vector  finite  elements;  (e)  the  mapping 
of  the  problem  to  a  matrix  equation  using  Galerkin’s  method;  (0  the  solution  of  the  system  of  equations 
using  LU  decomposition  or  the  conjugate  gradient  method;  and  finally,  (g)  the  implementation  of  all  of 
the  above  in  a  general-purpose  computer  code. 

There  are  numerous  possible  extensions  and  improvements  of  the  work  presented  here.  The  most 
important  of  these  are:  (a)  improving  the  efficiency  of  the  matrix  solution  by  using  preconditioning  in 
conjunction  with  the  conjugate  gradient  method;  and  (b)  implementation  of  a  "combined  source"  analogue 
for  the  aperture  field  continuity  to  prevent  the  kind  of  resonance  effects  observed  with  the  coax -rectangu¬ 
lar  launcher  problem.  A  few  others  that  would  be  quite  straightforward  code  modifications  are:  (a) 
increasing  the  number  of  ports;  and  (b)  allowing  ports  to  be  located  in  the  side  and/or  top  walls  of  the 
cavity.  Others  that  are  desirable,  but  will  require  some  theoretical  development  as  well,  are  anisotropic 
materials  and  active  devices.  However,  as  stated  in  the  introduction,  this  project’s  goal  is  to  apply  the 
hybrid  finite  element  method  to  phased  array  antenna  analysis,  and  our  work  for  the  immediate  future 
will  continue  in  that  direction. 
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Appendix  A 

The  Electric  Field  Functional 


Al.  Variational  Principle  vs  Weak  Form 

Publications  dealing  with  electromagnetic  finite  element  applications  usually  begin  with  one  of 
two  forms  of  functionals  '  with  little  or  no  discussion  as  to  why  one  is  used  and  not  the  other.  This 
appendix  discusses  their  origins  and  the  circumstances  in  which  each  is  appropriate.  The  conclusion  is 
that  the  differences  in  the  results  of  the  two  methods  are  inconsequential  for  typical  electromagnetic 
radiation  and  scattering  problems  even  though  only  one  of  them  is  considered  to  be  rigorously  justified. 

The  first  alternative  is 


F,{E) 


—  V'kE'  V  xE  -  kQe^E'  E 

Mr 


^  jj  (ExH)'  fids 

v 


(Al) 


(see,  for  example,  Jin  &  Volakis  (24])  where  fl  is  the  volume  region  over  which  the  unknown  field  is 


'  A  functional  is  an  integral  whose  integrand  includes  an  unknown  function.  It  maps  the  function 
into  a  number,  that  is,  the  result  of  the  integration  is  a  single  number,  as  opposed  to  an  integral  equation, 
which  maps  the  function  into  another  function. 
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to  be  found  and  F  is  its  enclosing  boundary.  The  subscript  v  denotes  that  this  functional  is  the  variational 
principle?  The  second  alternative  is 


(W^H)-Ads 


(A2) 


(see,  for  example,  D’Angelo  &  Mayergoyz  (7J).  Here,  W  is  a  trial  function  whose  form  is  yet  to  be 
determined.  The  subscript  w  on  F  denotes  the  fact  that  this  is  the  weak  form. 

The  variational  principle  F^,  has  been  deduced  from  a  general  from  of  energy  functional  |2S] 
and  is  used  in  conjunction  with  the  Rayleigh-Ritz  principle  to  form  a  system  of  equations.  The  variational 
principle  is  usually  avoided  in  situations  where  the  problem  is  not  self-adjoint  because  the  functional  F^, 
cannot  then  be  proven  to  be  stationary  about  the  solution  to  the  operator  equation.  The  weak  form  is 
derived  more  directly  by  simply  taking  the  inner  product  of  the  operator  equation  with  the  trial  function. 
It  will  be  used  in  conjuction  with  the  method  of  weighted  residuals  to  form  the  system  of  equations. 
When  Galerkin’s  method  (a  specialization  of  weighted  residuals)  is  used,  the  two  forms  may  give  exactly 
the  same  system  of  equations.  But,  to  use  Galerkin’s  method,  the  expansion  function  that  is  admissible 
in  the  original  problem  must  also  be  admissible  in  the  adjoint  problem. 

The  following  section  will  discuss  the  meaning  of  the  adjoint  problem  and  what  the  conditions 
are  for  self-adjointness  as  applied  to  the  Helmholtz  equation.  The  third  section  will  then  show  that  typical 
waveguide  continuity  conditions  represent  non-self-adjoint  boundary  conditions.  The  fourth  will  show 
that  under  the  assumption  that  Galerkin’s  method  is  applicable,  the  two  formulations  generate  identical 
systems  of  equations. 


^  A  variational  statement  may  be  either  a  variational  principle  or  a  weak  form. 
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A2.  The  Adjoint  Problem 


To  decide  whether  to  use  Eq.(Al)  or  Eq.(A2)  one  must  know  whether  or  not  the  operator 
equation  is  self-adjoint.  If  it  is  not,  then  the  properties  of  the  adjoint  problem  must  be  fouiid  in  order 
to  ensure  that  the  trial  functions  are  capable  of  representing  its  solutions. 

The  operator  equation  is  the  time-harmonic,  source-free  Helmholtz  equation  for  the  electric  field 
in  a  linear,  isotropic,  inhomogeneous  region: 


L(£)  =  Vx-Lvxf-Kt  £  =  0 

Its  inner  product  with  an  arbitrary  complex  function  VV  is 


<L(E),wy 


J.  Vx£- 


kQe 


W*  dv  =  0 


(A3) 


(A4) 


Using  a  Green’s  identity  (integration  by  parts)  twice,  we  shift  all  of  the  derivatives  ^'om  £  to  W: 


</-(£),  W> 


iii[i 


Vx£*  VxW  W* 


dv 


II  ±W*  XV  xE>  fids 
r 


(A5) 


From  the  definition 


Vx-LvxW-kle*  W 


dv 


Q 

-  II  ±[VV*  xVxE-ExVxW*]-  fids 


(A6) 


<£(£),  IV>  =  iE,L°{V/))  (A7) 

it  is  evident  that  the  term  in  brackets  in  the  first  integral  of  (A6)  is  £",  the  adjoint  operator.  It  is  simply 
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the  Helmholtz  equation  with  the  constitutive  parameters  replaced  by  their  complex  conjugates.  It  is  now 
clear  that  W  must  be  in  the  domain  of  the  adjoint  operator. 

The  definition  of  self-adjointness  is  L=L°,  which  obviously  cannot  be  true  if  the  problem  includes 
lossy  materials.  But  it  also  depends  on  whether  the  boundary  conditions  are  such  as  to  make  the  surface 
integral  in  Eq.(A6)  vanish,  that  is, 


J|  Ie^^*  xV  xE‘  nds  =  II  -Lfx  Vx£^*  • 
r  "  r 


(A8) 


^  ^E^*  xH^  Ms  -  II  £  X *  • 

'  >‘r  r 


fids 


(A9) 


Suppose  that  all  lossy  magnetic  materials  are  confined  inside  fl  so  that  along  the  boundary  F,  the  perme¬ 
ability  is  entirely  real.  Tlien  Eq.(A9)  is  satisfied  when  E^=E*  and  Recall  that  the  time-har¬ 

monic  and  time-dependent  fields  are  related  by  [Reference  26,  p.l5] 


E{x,y,z,i)  =  v^ReCEt’^"') 

(boldface  represents  the  time-dependent  quantity;  the  expression  for  magnetic  field  is  identical).  Evaluat¬ 
ing  Eq.(AlO)  at  any  point  t  along  F,  with  Eq  denoting  E{f)\ 


E(r,{)  =  \/^[^Re(£o)coswr  -  ImCfQjsinwr 


(All) 


and  assuming  that  B‘=E* 
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E“(r,r)  =  v^Re  £0*^^“']  =  y^[Re(£o)cosw/  +Im(£o)  sinu/j 

=  ^2  Re[£o<'‘^“']  (A12) 

=  ECr,-t) 

In  other  words,  the  adjoint  fields  are  time-reversed  versions  of  the  original  fields.  They  carry  power 
across  the  boundary  in  the  opposite  direction  and  they  encounter  materials  that  have  gain  instead  of  loss. 
This  is  the  physical  interpretation  of  the  adjoint  problem,  and  is  consistent  with  the  property  that  if  the 
operator  L  is  causal,  then  the  operator  £“  is  anti-causal  (Reference  27,  p.3561. 

Notice  that  if  the  boundary  F  is  comprised  entirely  of  perfect  electric  and  perfect  magnetic 
conductors,  then  Eq.(A9)  is  always  satisfied  because  there  cannot  be  any  transfer  of  power  across  such 
boundaries.  Thus  we  are  led  to  suspect  that  the  boundary  conditions  that  will  cause  non-self-adjointness 
are  those  open  boundaries  where  conditions  of  field  continuity  are  to  be  enforced.  The  next  section  will 
demonstrate  that  this  is  true  for  the  waveguide/cavity  apertures  that  are  considered  in  the  main  body  of 
this  report. 

A3.  Continuity  Conditions  for  Waveguide  Apertures 

For  the  waveguide  continuity  problem,  the  boundary  F  reduces  to  the  waveguide  aperture  F^^. 
Taking  ^=E*  and  the  following  form  is  equivalent  to  Eq.(A9): 

II  £*  -(flxWjdj  =  II  £.(/}xW*)(/5 

(again  assuming  that  Ur 's  real  along  F).  Consider  a  waveguide  whose  axis  is  the  z  axis  and  which  joins 
the  volume  fl  through  an  aperture  in  its  end  wall  at  z=0.  It  is  assumed  to  be  match-terminated  at  z<<0. 
The  field  in  the  aperture  due  to  the  dominant  mode  and  its  complex  conjugate  are 
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(A  14) 


E  =  «o(‘"Co) 

E*  =goO*Co) 


H  =  f^«o^o(l-Co) 
H*  ^  t^goYoi^-Co) 


(A15) 


where  Ig  ^  transverse  mode  function,  Yq  is  the  modal  admittance  and  Cq  is  an  unknown  coefficient. 
The  outward  normal  to  0  is  =  so  the  left  and  right  sides  of  Eq.(A13)  are 


E 


•(AxH)ds  =  Yo{l*Co){l- 


~  ds 


(A  16) 


II  E-{flxH*)ds  =  Tod^CoXl-Co  )||  (A17) 

respectively.  They  are  only  equal  if  Cg  is  real,  which  is  generally  false.  Therefore,  the  continuity  condi¬ 
tion  across  a  waveguide  aperture  will  render  the  operator  equation  non  self-adjoint;  the  functional  Eq.(A  1) 
is  not  appropriate  even  when  there  are  no  lossy  materials;  and  the  weak  form  should  be  used. 

A4.  Galerkin’s  Method  vs  Rayleigh-Ritz 

In  the  main  body  of  this  report,  the  finite  element  method  is  used  to  produce  a  system  of  equa¬ 
tions  from  the  weak  form  functional.  The  field  is  represented  by  a  summation  of  unknown  complex 
coefficients  with  known,  linear  vector  functions: 

£  =  £  =  52  e^^,{x,y,z)  (A18) 

For  the  expansion  functions  to  be  admissible,  they  must  be  in  the  domain  of  the  functional.  Linear 
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functions  meet  that  requirement  since  their  first  derivatives  are  continuous  (integrable). 

The  residual  error  is  defined  as  R=L(£).  N  weighting  functions  will  be  chosen,  each  of  which 
is  orthogonal  to  the  residual  so  that  {R,w^)=Q,  giving  the  system  of  N  equations 


N 

0  ■ 

5=1 


I  _  'y  _ 

—  V  X  •  V  X 
Mr 


dv 


(A19) 


*  J^oVo  \\  >  '■=1,2 . N 


The  choice  satisfies  the  orthogonality  requirement,  but  it  is  also  necessary  that  be  an  admissible 

expansion  function  for  the  adjoint  electric  field.  Starting  with  the  adjoint  operator  equation,  forming 
(£■,/<“(  W))  and  using  the  Green’s  identity  once  gives 


<£,L“(Vf)> 


±  V  X  W* 
Mr 


VxE-kle^W* 


E  dv 


(( 

f* 


J-ExVxW* 

Mr 


•  fids 


(A20) 


which  makes  it  evident  that  expansion  functions  for  W must  have  continuous  first  derivatives.  Therefore, 
the  same  expansion  functions  are  admissible  in  both  the  original  and  the  adjoint  problems.  It  is  this  fact 
that  allows  Galerkin’s  method  to  be  employed. 

Consider  the  variational  functional,  Eq.(Al),  with  the  expansion  Eq.(A18)  for  the  electric  field; 


N  N 

r=l  j=l 

N 


0  L 


Mr 


dv 


(A21) 


^rjj  >'=1,2 . 

r=l  r 


N 


The  Rayleigh-Ritz  method  equates  the  stationary  point  of  F„  to  the  minimization  of  the  above  with  respect 
to  each  of  the  coefficients  e^,  that  is. 
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5F,(E)  =  0  « 


(A22) 


dF,(£) 

~de^ 


=  0, 


Vr 


Carrying  out  the  partial  derivatives  in  Eq.(A21)  gives 


N 

0  =  E 

0  L 


J_  V  X  •  V  X  -  itg 
Mr 


(A23) 


Jl^oVo  I  r=l,2,...,N 


which  is  identical  to  Eq.(A  19)  with  replaced  by  This  shows  that  the  systems  of  equations  resulting 
from  the  variational  principle  and  from  the  weak  form  are  identical.  Thus,  under  at  least  some  circum¬ 
stances  the  distinction  between  the  two  forms  is  inconsequential,  the  adjoint  form  may  be  ignored,  and 
the  variational  principle  may  be  used  in  situations  where  it  is  not  rigorously  Justified. 
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Appendix  B 

Coaxial  and  Circular  Waveguide  Mode  Functions 


B.l.  COAXIAL 

In  most  instances  the  properties  of  circular  coaxial  waveguide  are  governed  by  the  TEM  mode 
because  higher-order  modes  are  strongly  attenuated.  Near  discontinuities,  however,  those  higher  order 
TE  and  TM  modes  contribute  significantly  to  the  fields,  and  so  their  effects  must  be  included  to  accurate¬ 
ly  predict  the  properties  of  discontinuities.  This  appendix  discusses  the  nature  of  coaxial  modes  and  gives 
formulas  defining  the  orthonormal  mode  functions  (extracted  or  derived  from  Marcuvitz  \2\).  It  also 
describes  a  set  of  computer  subroutines  that  implement  the  formulas,  as  well  as  some  calculations  that 
validate  their  accuracy. 

Bl.l.  Mode  Functions 

The  circular  coaxial  waveguide  has  inner  and  outer  radii  a  and  b,  respectively.  The  formulas  that 
follow  will  be  given  in  terms  of  polar  coordinates,  with  the  waveguide  axis  coinciding  with  the  z  axis. 
The  higher  order  modes  have  indices  m  and  n,  which  are,  respectively,  the  number  of  cycles  in  the 
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circumferential  direction,  and  the  number  of  half-cycles  in  the  radial  direction.  For  each  combination 
m.n  there  are  four  distinct  mode  types:  TE  or  TM;  and  cos(/n<^)  or  sin(/n<^).  All  four  must  be  included 
to  represent  an  arbitrary  field.  The  mode  types  will  be  represented  in  the  following  by  the  subscripts/? 
and  q  where  /?=  1  for  TE  or  2  for  TM,  and  q=  1  for  cos(m</»)  or  2  for  sin(m4>). 

Waveguide  fields  are  represented  in  terms  of  transverse  (no  t  component)  mode  functions 
The  transverse  component  of  an  arbitrary  field  with  propagation  in  the  +z  direction  is  represented  as  a 
sum  of  modes  with  unknown  excitation  coefficients 

>  E  E  E  E  (BD 

p=l  </=l  m=0  n=rt| 

The  starting  index  for  n  is  nj  =  2  for  p=  1  and  m=Q  or  nj  =  l  otherwise,  because  there  is  no  TEqj  mode  - 
the  first  TE  mode  is  mn=02.  The  propagation  constants  are 
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,  k  <  (evanescent) 

where  the  cutoff  wavenumber  k^  is  the  /t’th  root  of  one  of  the  characteristic  equations 

{TE) 


{TM) 

The  higher-order  mode  functions  are  given  by  the  following  ten  equations; 
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g2lmn  =  -P^mn22„„(*„„p)COS(/«<i.)  +  $  ^  Z2„„(k„„p)  sin(m<i>) 
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The  TEM  mode  function  is 

«»■? 

Finally,  it  is  convenient  to  represent  the  transverse  magnetic  fields  as 


(BIS) 
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where  is  a  modal  admittance 
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(B17) 


The  dominant  mode’s  admittance  is  simply  YQ=Hri.  A  subtle  but  important  point  is  that  this  modal 
admittance  is  not  the  same  as  the  characteristic  admittance  |28).  (The  characteristic  admittance  of  the 
dominant  mode  in  a  coaxial  cable  is  27r/[7;ln(^/a)|.) 


B1.2.  Computations 

The  IMSL®  library  routines  are  used  to  compute  Bessel  functions.  Additional  routines  were 
created  to  implement  formulas  for  the  Bessel  function  derivatives. 

Figures  B1  and  B2  show  the  TE  and  TM  characteristic  equations  for  various  orders  m  (the  ratio 
ft/a=7/3  is  typical  of  50  fl  characteristic  impedance).  The  problem  of  finding  the  roots,  indicated  by  the 
circles,  turned  out  to  be  non-trivial.  Abramowitz  &  Stegun  (Reference  10,  p.3741  give  three  terms  of 
a  series  for  the  roots.  Evidently,  more  terms  are  needed  for  orders  higher  than  m  =  3  and/or  n  =  4,  because 
in  those  instances  the  calculated  root  was  not  close  enough  to  the  actual  root  for  a  Newton’s  method 
iteration  to  converge  to  the  actual  root.  Unfortunately  it  is  not  possible  to  deduce  what  the  remaining 
terms  in  the  series  are  from  the  ones  that  are  given.  It  is  probably  no  coincidence  that  the  tables  of  roots 
given  in  Marcuvitz  [Reference  2,  pp. 74,78]  contain  blanks  for  those  values  where  the  3-term  series  fails. 

Marcuvitz  gives  approximate  formulas  for  nt=0  that  are  fairly  accurate  for  all  n.  The  procedure 
that  was  adopted  uses  this  as  a  seed  for  Newton’s  method  in  order  to  find  the  0,n  root.  That  is  used  in 
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turn  as  a  seed  to  find  the  1  ,n  root;  and  so  on  up  to  the  desired  m,n  root.  This  method  has  reliably  found 
all  of  the  roots  up  to  m.n=6,l0.  Figures  B3  and  B4  are  results  obtained  using  that  routine  to  find  the 
cutoff  frequencies  of  higher-order  modes  in  APC-7  coax  (fl=  1.502  nun,  i>=3.5  mm).  As  expected,  the 
first  mode  to  become  propagating  is  the  TEj,  mode,  but  its  cutoff  frequency  of  19.5  GHz  is  well  above 
the  specified  operating  limit  of  18  GHz  for  the  APC-7  coaxial  connectors. 

Once  the  cutoff  wavenumbers  are  known  it  is  a  tedious  but  straightforward  matter  to  find  the 
mode  functions  using  Eqs.(B5)-(B14).  Figures  B5  and  B6  are  the  results  of  example  calculations.  They 
are  contour  plots  of  the  electric  field  lines  for  two  higher  order  modes.  These  results  agree  with  figures 
in  Marcuvitz  [Reference  2,  pp. 76,79). 

The  important  feature  of  orthonormality  was  also  verified  by  numerically  integrating  function 
products  for  various  values  of  the  indices  p,q,m,n.  The  integral  of  any  pair  of  mode  functions  over  the 
entire  waveguide  cross  section  should  equal  1  if  the  two  functions  are  the  same,  and  should  equal  0  other¬ 
wise.  The  success  of  this  test  for  all  of  the  30  mode  pairs  attempted  verifies  the  accuracy  of  the  normal¬ 
ization  factors  (B13)  and  (B14). 

B2.  CIRCULAR  WAVEGUIDE 

The  procedures  used  for  circular  waveguide  functions  are  substantially  the  same  as  those  used  for 
coaxial  waveguides.  There  are  again  four  types  of  mode  functions,  TE  and  TM  with  sin(m0)  or  cos(m<f») 
dependence.  The  cutoff  wavenumbers  for  the  m,n  modes,  are  the  n’th  nonzero  roots  of  J^(ka)  for 
TE  {p~  1),  or  of  J„{ka)  for  TM  (p=2),  where  a  is  the  waveguide’s  inside  wall  radius.  The  same  root- 
finding  procedure  discussed  in  the  preceding  section  was  used  to  determine  the  cutoff  wavelengths.  The 
mode  functions  are  as  follows: 
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Figure  B3.  APC-7  Mode  Cutoff  Frequency  vs.  Index  n 
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Figure  B4.  APC-7  Mode  Cutoff  Frequency  vs.  Index  m 
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Appendix  C 

Rectangular  Waveguide  Mode  Functions 
and  Inner  Product  Integrations 

The  rectangular  waveguide  mode  functions  aredefined  for  a  waveguide  whose  side  walls  are  parallel 
to  the  X  and  y  axes,  with  lengths  a  and  b  respectively.  The  inside  cross  sectional  area  is  confined  to  the 
first  quadrant  so  that  two  of  the  inner  walls  are  in  thex-z  andy-z  planes.  The  matrix  terms  relating  interior 
electric  field  to  waveguide  modes  are  given  in  terms  of  products  of  integrals  of  the  form 

^  rmnp  ~  'Smnp^^  (C.l) 

□ 

where  □  denotes  the  open  portion  of  the  waveguide  end  wall  aperture,  's  the  vector  expansion  function 
associated  with  edge  r  and  is  the  transverse  mode  function  for  the  mn  mode  with  p=\  for  TE  and 
p=2  for  TM.  Since  is  only  defined  over  the  triangles  adjacent  to  edge  r  (2  at  most),  the  region  of 
integration  reduces  to  those  triangles.  Over  asingletrianglewith  index  suppose  that  the  edge  r’sendpoints 
are  the  local  nodes  1  and  2.  Then  substituting  the  definition  of  the  expansion  function  Eq.(20)  gives 
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*!lp  -  (  (/i  ^f2  -h  -is 


(C2) 


where /|  and  fn  are  the  linear  scalar  expansion  functions  associated  with  nodes  I  and  2,  respectively. 
Transforming  to  homogeneous  coordinates  gives /|  =  and^  =  /t.  dx  dy  =  2A  dt^  rf/o 


’/l  -  5!;(<7',2*)’r„)  (C3) 


where  A  is  the  triangle  area.  The  matrix  T  gives  the  coordinate  transformation  as  in  Eq.(l6).  The  inverse 
of  this  transform  gives  x  and  y  in  terms  of  /,  and  /2  as 


Jt  =  aQ  +  a,  r,  +  02  ^2 

(C5) 

y  =  00* 0^t^  *02h 

(C6) 

The  transverse  mode  functions  and  normalization  coefficients  are  given  by  Marcuvitz  (Reference 
15,  pp.57-601; 
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(Cll) 


■  -1/2 

n~a  ^  m~  b  (C12) 

4b  4  a 

where  is  Neumann’s  number. 

Let  the  coefficients  and  jj-  represent  the  following  combinations  of  a-  and  from  Eqs.(C5) 
and  (C6); 
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which  allows  the  trigonometric  functions  in  Eqs.(C7)  and  (C8)  to  be  expressed  in  the  convenient  form 


cos(^)sinj^j  =  isin(^o-*-{,r,  ■♦•^2^2)"isin(T/ot-Tj,/, +1/2^2) 


(CIS) 


i(i!l^)cos(^)  =  :Jsin(^o  +  ^,r,+^2^2)'":7Sin(7;o  +  7J,/, +7/2^2) 


(Cl  6) 


Let  f/j  denote  the  integral 


//j  =  I|  t^dl^  j  sin(^o  +  i,  r, +i2^2)‘^^2 

~  0  0 


(Cl  7) 


//^  is  identical  with  tj-  replacing  everywhere.  and  are  the  same,  but  with  and  ^2  or  t/i  and 
7/2  reversed.  Then  finally,  the  inner  product  integral  may  be  written  as 
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(Cl  8) 


-  T,,C„,n^'<r>*'v>  -  ^i3C,,„(«i'*H,')| 

(The  transformation  of  the  integral  to  homogeneous  coordinates  introduces  a  scale  factor  of  2/4  that  cancels 
the  1/2/4  in  Eqs.(C3)  and  (C4).)  The  integral  in  Eq.(C17)  may  be  carried  out  in  closed  form,  but  it  must 
be  evaluated  separately  for  each  of  the  following  cases  (a)  through  (e): 

(a)  general  case: 


H, 


—  [  cos(^o  i  1 )  -  cos^o  +  sin(4o  f  i )] 
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(b)  ^2  5*^03  i  1 I  ^2'^  1 1  '  (((’<■  small  values  of  {2'i  i>  Eq.(C.  19a)  is  numerically  unstable, 

susceptible  to  large  roundoff  error  and  overflow); 
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(c)  U2l^C 


//|  =  — !-t  |(4^,  ■^6^2"^r^2)cos(^o) -(4^, +6^o)COS(^o-^^l) 
-(2^f  +4^,^2)sin(^o)  -2({f  +  ^,^2)sin(^o-^^i)} 

(d)  Uil<l: 
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(e) 


I  sin^o‘"(y  ■"-^)cos(^o)  I 

Note  that  the  last  four  are  correct  in  the  limits  as  those  quantities  specified  as  much  less  than  unity  go 
to  zero.  The  computer  code  implementation  gives  precedence  to  the  five  formulas  in  reverse  order,  checking 
for  condition  (e)  first,  and  executing  (a)  only  when  none  of  the  others  (d),  (c)  or  (b)  are  true. 
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MISSION 

OF 

ROME  LABORATORY 

Rome  L,ahoratory  plans  and  executes  an  interdisciplinary  program  in  re¬ 
search,  development,  test,  and  technology  transition  in  support  of  Air 
Force  Command,  Control,  Communications  and  Intelligence  (C^I)  activities 
for  all  Air  Force  platforms.  It  al^  executes  selected  acquisition  programs 
in  several  areas  of  expertise.  Technical  and  engineering  support  within 
areas  of  competence  is  provided  to  BSD  Program  Offices  (POs)  and  other 
BSD  elements  to  perform  effective  acquisition  of  C^I  systems,  h  additicm, 
Rome  Laboratory's  technology  supports  other  AFSC  Product  Divisions,  the 
Air  Force  user  community,  and  other  DOD  and  non-DOD  agencies.  Rome 
Laboratory  maintains  technical  competence  and  research  programs  in  areas 
including,  but  not  limited  to,  communications,  command  and  control,  battle 
management,  intelligence  information  processing,  computational  sciences 
and  software  producibility,  wide  area  surveillance/sensors,  signal  proces¬ 
sing,  solid  state  sciences,  photonics,  electromagnetic  technology,  super¬ 
conductivity,  and  electronic  reliability/maintainability  and  testability. 


